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ABSTRACT 


This  technical  report  summarizes  the  image  processing  re¬ 
search  activities  performed  by  the  University  of  Southern  California 
during  the  period  of  1  March  1974  to  31  August  1974  under  Contract 
No.  F08606-72-C-0008  with  the  Advanced  Research  Projects  Agency, 
Information  Processing  Techniques  Office. 

The  research  program,  entitled,  "Image  Processing  Research," 
has  as  its  primary  purpose  the  analysis  and  development  of  techniques 
and  systems  for  efficiently  generating,  processing,  transmitting,  and 
displaying  visual  images  and  two  dimensional  data  arrays.  Research 
is  oriented  toward  digital  processing  and  transmission  systems.  Five 
task  areas  are  reported  on:  (1)  Image  Coding  Projects,  the  investigation 
of  digital  bandwidth  reduction  coding  methods;  (2)  Image  Restoration  and 
Enhancement  Projects:  the  improvement  of  image  fidelity  and  presentation 
format;  (3)  Image  Data  Extraction  Projects:  the  recognition  of  objects 
within  pictures  and  quantitative  measurement  of  image  features;  (4)  Image 
Analysis  Projects,  the  development  of  quantitative  measures  of  image 
quality  and  analytic  representation;  (5)  Image  Processing  Support  Projects. 
development  of  image  processing  hardware  and  software  support  systems. 


I 


PROJECT  PARTICIPANTS 


Project  Director 
William  K.  Pratt 

Research  Staff 
Harry  C.  Andrews 
Lee  D.  Davisson 
Werner  F rei 
Ali  Habibi 
Ernest  L.  Hall 
Ronald  S.  Hershel 
Anil  K.  Jain 
Richard  P.  Kruger 
Nasser  E.  Nahi 
Guner  Robinson 
Alexander  A.  Sawchuk 
Lloyd  R.  Welch 

Support  Staff 
Angus  B.  Cossey 
Toyone  Mayeda 
Paula  Meeve 
James  M.  Pepin 
Easter  D.  Russell 
Mark  A.  Sanders 
Ray  Schmidt 
Joyce  Seguy 
Dennis  Smith 
George  Soen 
Wai  Szeto 

Florence  B.  Tebbets 


Students 
Ben  Britt 
Steve  Dashiell 
Faramarz  Davarian 
Roy  M.  Glantz 
Henry  Glazner 
Michael  Huhns 
Steve  Hui 

Mohammad  Janshahi 
Mohser.  Khalil 
San  Kwok 
Alan  Larson 
Paul  Liles 
Robert  Liles 
Eduardo  Lopez 
Simon  Lopez-Mora 
Clanton  Mancill 
Nelson  Mascarenhas 
David  Miller 
Firouz  Naderi 
Michael  Patton 
Mohammad  Peyrovian 
Stewart  Robinson 
Robert  Wallis 
Pamela  Welch 


■li- 


/ 


TABLE  OF  CONTENTS 


!•  Research  Project  Overview  j 

2.  Research  Project  Activities  2 

3.  Image  Coding  Projects  2 

3.1  Adaptive  Dual  Mode  DPCM/Deltamodulation  3 

Image  Coding  Techniques 

3.2  Orthogonal  Transform  Coding  of  Moving  Pictures  7 

3.  3  Quantization  Error  Reduction  for  Image  Coding  16 

4.  Image  Restoration  and  Enhancement  Projects  27 

4.1  Restoration  for  Space  Variant  Aberrations  28 

4.2  Image  Restoration  in  the  Eigenspace  of  the  Degradation  32 

4.  3  Pseudoinverse  Method  of  Bounded  Image  Restoration  34 

4.4  A  Fast  Pseudoinverse  Image  Restoration  36 

4.5  Spline  Function  Image  Restoration  45 

4.6  Nonuniform  Sampling  of  Observation  Space  48 

4.7  Histogram  Exponentiation 

5.  Image  Data  Extraction  Projects  59 

5.1  Fourier- Bessel  Method  for  Transverse-Axial  59 

5.2  Nonlinear  Optical  Image  Processing  with  Halftone  Screens  65 

6.  Image  Analysis  Projects 

6.1  A  Quantitative  Model  of  Color  Vision  69 

6.  2  A  New  Look  at  the  Visual  System  Model  84 

7.  Image  Processing  Hardware  and  Software  Projects  102 

7.1  Hardware  Projects  102 

7.2  Software  Projects  jq3 

8.  Image  Processing  Facilities  IO5 

8.1  Image  Processing  Laboratory  105 

8.2  Engineering  Computer  Laboratory  113 

8.  3  Optical  Information  Processing  Laboratory  117 

8.4  Biomedical  Image  Processing  Laboratory  119 

9.  Publications  12q 

-iii- 


This  report  describes  the  progress  and  results  of  the  University 
of  Southern  California  image  processing  research  study  for  the  period  of 
1  March  1974  to  31  August  1974.  The  image  processing  research  study 
has  been  subdivided  into  five  projects: 

Image  Coding  Projects 

Image  Restoration  and  Enhancement  Projects 
Image  Data  Extraction  Projects 
Image  Analysis  Projects 
Image  Processing  Support  Projects 

In  image  coding  the  orientation  of  the  research  is  toward  the  development 
of  digital  image  coding  systems  that  represent  monochrome  and  color  im¬ 
ages  with  a  minimal  number  of  code  bits.  Image  restoration  is  the  task 
of  improving  the  fidelity  of  an  image  in  the  sense  of  compensating  for  im¬ 
age  degradations.  In  image  enhancement,  picture  manipulation  processes 
are  performed  to  provide  a  more  subjectively  pleasing  image  or  to  convert 
the  image  to  a  form  more  amenable  to  human  or  machine  analysis.  The 
objectives  of  the  image  data  extraction  projects  are  the  registration  of  i.n- 
ages,  detection  of  objects  within  pictures  and  measurements  of  image  fea¬ 
tures.  The  image  analysis  projects  comprise  the  background  research  ef¬ 
fort  into  the  basic  structure  of  images  in  order  to  develop  meaningful  quan¬ 
titative  characterizations  of  an  image.  Finally,  the  image  support  projects 
include  research  on  image  processing  computer  languages  and  the  develop¬ 
ment  of  experimental  equipment  for  the  sensing,  processing,  and  display  of 
images. 

The  next  section  of  this  report  summarizes  some  of  the  iesearch 
project  activities  during  the  past  six  months.  Sections  3  to  7  describe  the 
research  effort  on  the  projects  listed  above  during  the  reporting  period. 

A  capsule  description  of  the  physical  facilities  of  the  USC  Image  Proces¬ 
sing  is  contained  in  Section  8.  Section  9  is  a  list  of  publications  by  project 
members. 


2.  Research  Project  Activities 


Significant  research  project  activities  of  the  past  six  months  are 
summarized  below: 

Summer  Short  Courses.  One  of  the  vehicles  for  the  transfer  of  ARPA  spon¬ 
sored  image  processing  research  technology  to  the  Federal  and  industrial 
communities  has  been  through  intensive  one  and  two  week  short  courses. 

For  the  past  four  years  the  University  has  offered  a  one  week  Summer 
short  course  in  Mathematical  Pattern  Recognition  and  a  two  week  course 
in  Digital  Image  Processing.  In  addition,  this  Summer,  a  third  short 
course  of  one  week  duration  was  initiated  on  Optical  Processing.  The 
alumni  of  these  courses  now  numbers  in  the  hundreds. 

Picture  Coding  Symposium.  Professor's  William  K.  Pratt,  Ali  Habibi, 
and  Werner  F rei  attended  the  Picture  Coding  Symposium  held  in  Goslar, 
West  Germany  on  26-28  August,  1974.  Approximately  150  persons  from 
throughout  the  free  world  attended  the  meeting  which  was  devoted  to  ses¬ 
sions  on:  human  observer  characteristics,  intraframe  image  coding,  inter¬ 
frame  image  coding,  color  image  coding,  and  multi- spectral  data  coding. 

At  the  conference  a  contest  was  held  to  choose  the  best  algorithms  for  pic¬ 
ture  coding  developed  during  the  past  18  months.  Dr.  Ali  Habibi  of  USC 
was  awarded  a  prize  for  the  best  coding  algorithm  for  monochrome  image 
coding  with  1.  0  bits /pixel. 

Special  Issues.  During  the  past  six  months  two  special  issues  relating  to 
image  processing  were  published.  Professor  Harry  C.  Andrews  of  USC 
was  editor  of  the  IEEE  Computer  journal  May  1974  special  issue  on  "Com¬ 
puter  Image  Processing."  The  Society  of  Photo-Optical  Instrumentation 
Engineers  devoted  its  May  1974  issue  of  Optical  Engineering  to  optical  and 
digital  image  processing  under  the  editorship  of  Professor  Alexander  A. 
Sawchuk  of  USC. 
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3.  Image  Coding  Projects 

The  effort  in  image  coding  is  directed  toward  the  research  and 
development  of  image  coding  systems  providing  a  transmission  bit  rate 
reduction  and  tolerance  to  channel  errors.  Coding  systems  are  under 
investigation  for:  monochrome  and  color  imagery:  slow  scan  and  real 
time  television:  and  information  preserving  and  controlled  fidelity 
operation.  Results  of  this  research  study  during  the  past  six  months 
are  summarized  here  and  presented  in  detail  in  subsequent  sections. 

The  first  report  concerns  a  performance  study  of  several 
adaptive  linear  predictive  image  coding  systems.  One  such  system 
which  adaptively  codes  by  DPCM  or  deltamodulation  dependent  upon 
picture  quality  has  been  found  to  provide  high  fidelity  coding  down  to 
about  2.0  bits/pixel. 

The  next  report  describes  an  application  of  Hadamard  and  Slant 
transforms  to  interframe  image  coding.  Good  results  are  obtained  at 
an  average  of  much  less  than  one  bit/pixel. 

In  the  final  report  simulation  results  are  given  tor  transform 
coding  spectrum  interpolation  in  which  receiver  post-processing  is 
employed  to  reduce  quantization  error  effects  on  transform  coefficients. 
Reduction  in  mean  square  error  of  about  30%  can  be  obtained  by  this 
technique. 

3.  1  Adaptive  Dual  Mode  DPCM/Deltamodulation  Image  Coding 
Techniques 

William  K.  Pratt 

Standard  one  bit  per  sample  deltamodulation  provides  a  consider¬ 
able  bandwidth  reduction,  but  suffers  from  the  complementary  problems 
of  image  granularity  in  smooth  areas  of  an  image  and  slope  overload  in 
image  regions  of  rapid  brightness  change.  These  problems  may  be 
alleviated  by  quantizing  the  difference  signal  with  more  levels  as  in  3  bit 
per  sample  DPCM,  but  of  course,  the  bandwidth  reduction  is  sacrificed. 
A  compromise  technique  has  been  investigated  in  which  the  image  coder 
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operates  in  a  dual  mode:  DPCM  in  regions  of  high  image  activity  and 
deltamodulation  in  regions  of  low  activity.  The  switch  between  modes 
is  performed  adaptively  based  upon  tne  relative  image  activity. 

In  a  basic  deltamodulation  image  coding  system  a  prediction  of 
the  next  pixel  to  be  scanned  is  made  based  upon  the  previously  scanned 
pixel.  The  difference  between  the  actual  pixel  value  and  its  estimate  is 
quantized  to  two  levels  and  transmitted  as  a  binary  pulse  (one  bit  pixel 
code).  Figures  la  to  lc  contain  computer  simulation  photographs  of 
deltamodulation  coded  pictures  for  three  quantization  levels  'q  =  5%,  10%, 
20%)  and  a  previous  pixel  weighting  of  90%.  The  tradeoff  between  slope 
overload  error  and  granularity  error  is  readily  apparent  from  the  photo¬ 
graphs. 

With  DPCM  encoding  the  difference  signal  between  a  scanned 
pixel  and  its  previous  pixel  estimate  is  quantized  to  8  levels  '3  bit/pixel 
code)  which  are  nonlinearly  spaced  to  minimize  quantization  error.  Fig  ir? 

Id  contains  an  example  of  DPCM  coding  when  the  first  threshold  level  is 
3et  at  2.  5%  full  scale  and  a  90%  previous  pixel  weighting  is  employed. 

The  resultant  image  quality  is  satisfactory  for  most  applications. 

One  possibility  for  reducing  the  bit  rate  requirement  of  DPCM  is 
to  employ  a  dual  mode  coder  which  switches  from  DPCM  to  deltamodulation 
in  regions  of  nearly  constant  grey  level  and  from  deltamodulation  to  DPCM 
in  edge  regions.  A  simple  algorithm  for  the  switch  is  as  follows: 

(a.)  Switch  from  deltamodulation  to  DPCM  if  three  sequential  delta 
bits  are  of  the  same  sign 

(b. )  Switch  from  DPCM  to  deltamodulation  if  DPCM  quantizer  shifts 
from  smallest  positive  to  smallest  negative  quantization  level 
or  vice  versa. 

By  adjusting  the  quantization  levels  for  the  DPCM  and  deltamodulation 
quantizers  it  is  possible  to  control  the  relative  time  division  between 

states.  One  of  the  major  advantages  of  the  encoder  is  that  the  decoding 
may  be  performed  at  the  receiver  from  the  transmitted  code  without  any 
explicit  code  bits  required  to  designate  the  mode.  Figure  le  contains 
a  coded  image  using  tnis  algorithm  coded  at  about  2.0  bits/pixel.  In 
regions  of  little  image  activity,  image  quality  is  good,  but  mode  transition 
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(d)  DPCM 

Q  =2.  5%,  3  bits/pixel 


(a)  deltamodulation 
q  =  5%,  1  bit/pixel 


(e)  dual  mode 
2  bits /pixel 


(b)  deltamodulation 
q  =10%,  1  bit/pixel 


(c)  deltamodulation 
q  =  20%,  1  bit /pixel 


(f)  oversampled  dual/mode 
2  bits /pixel 


Figure  3.1-1.  Examples  of  deltamodulation,  DPCM 
and  dual  mode  coded  images 


errors  are  apparent  near  edges. 

Frei,  Schindler,  and  Velliger  [}]  have  suggested  a  dual  mode 
system  in  which  the  original  image  is  oversampled  by  a  factor  of  three 
to  permit  a  rapid  detection  of  the  delta  to  DPCM  mod  e  change.  The 
coding  logic  is  as  follows: 

(a.  )  Switch  from  delatamodulation  to  DPCM: 


After  three  sequential  delta  "ones"  insert  a  "zero"  marker  bit  and 
after  three  sequential  delta  "zeros"  insert  a  "one"  marker  bit. 

(b.)  Switch  from  DPCM  to  deltamodulation  if  DPCM  quantizer  shifts 

•ft 

from  smallest  positive  to  smallest  negative  quantization  level  or 
vice  versa. 

Figure  1  f  contains  a  coded  image  for  the  over  sampled  dual  mode 
coding  system  for  coding  at  two  bits  per  pixel.  The  oversampled  dual 
mode  system  is  superior  to  the  simple  dual  mode  system  described 
previously  in  terms  of  picture  quality,  however,  its  implementation 
requirements  are  greater. 

In  summary,  the  oversampled  dual  mode  DPCM/deltamodulation 
image  coding  system  has  proven  to  provide  good  quality  coded  images  at 
2.0  bits/pixel.  The  di sadvantage  of  the  coder  is  the  additional  complex¬ 
ity  as  compared  to  conventional  3  bit /sample  DPCM. 
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3.2  Orthogonal  Transform  Coding  of  Moving  Pictures 
Clifford  Reader 

The  efficiency  of  coding  a  sequence  of  moving  pictures  is  improved 
by  techniques  of  interframe  coding  which  partially  lemove  redundancy 
between  frames.  A  successful  type  of  coding  is  that  of  conditional  update 
operating  upon  the  differences  between  successive  frames.  Orthogonal 
transform  techniques  have  been  applied  to  the  conditional  update  process 
with  additional  advantages  over  conventional  systems. 

Although  interframe  codir.^  may  be  very  successful  in  reducing 
the  transmission  requirements  for  moving  picture  signals,  the  coder  is 
usually  complex  and  costly,  involving  in  expensive  frame  memory.  The 
object  of  the  orthogonal  transform  conditional  update  coder  is  to  minimize 
terminal  costs.  The  coder  performs  a  two  stage  process.  The  first 
stage  is  a  conventional  intraframe  coder,  removing  spatial  redundancy 
from  single  frames.  This  orthogonal  transformation  and  block  quan¬ 
tization  is  performed  over  sub-blocks  within  the  frame.  The  second, 
interframe  coding  stage,  then  applies  conditional  update  techniques  to 
the  differences  between  successive  transformed  and  quantized  sub-blocks. 
This  is  simply  achieved  by  setting  a  threshold  and  updating  a  sub-block  if 
the  energy  of  the  difference  signal  foi  that  sub-block  exceeds  the  thres¬ 
hold.  A  block  diagram  of  the  coder  is  shown  in  figure  1.  It  is  to  be 
noted  that  the  updating  process  takes  place  with  the  data  from  the  new 
frame  and  not  the  difference  information.  This  is  done  to  avoid 
problems  with  quantizing  the  difference  signal  which  is  derived  from 
the  nonlinear  intraframe  quantizer.  An  advantage  is  that  the  memory  is 
refreshed  with  new  data  and  not  an  estimate  of  the  new  data.  This  reduces 
the  problem  of  multiplicative  error.  The  principle  advantage  of  the 


-7- 


interframe  coder  is  derived  from  the  use  of  intraframe  coded  data. 

This  data  is  thus  presented  at  reduced  rate  and  the  memory  requirement 
of  the  coder  is  reduced.  If  the  interframe  coder  produces  a  compres¬ 
sion  of  M:1  then  the  size  of  the  deJays  must  be  1/256M  frames  (oius  a 
small  amount  to  allow  for  the  difference  signal  energy  thresholding  to 
take  place)  and  the  memory  size  will  be  1/M  -  1/256M  =  255/256M 
Frames.  The  overall  memory  (storage)  requirement  for  the  coder  is 
thus  reduced  by  almost  a  factor  of  M  compared  with  the  requirement  for 
conventional  coders. 

Quantitative  analysis  of  the  interframe  error  is  hampered  by  the 
nonlinearity  of  the  updating  process.  Qualitatively,  several  observa¬ 
tions  may  be  made.  In  the  absence  of  noise,  a  threshold  of  zero  would 
result  in  the  transmission  of  all  sub-blocks  posessing  a  finite  difference 
signal  energy  and  the  interframe  error  would  be  zero.  If  the  thres¬ 
hold  were  then  raised,  those  sub-blocks  which  contained  the  least 
amount  of  motion  would  not  be  updated.  Since  the  difference  signal 
energy  is  closely  related  to  the  mean  squared  error  made  when  a  sub¬ 
block  is  not  updated,  it  is  reasonable  to  assume  that  the  interframe  error 
should  be  proportional  to  some  function  of  the  update  threshold.  In  the 
presence  of  noise,  this  effect  will  be  modified.  The  interframe  error 
was  examined  for  five  pairs  of  frames  containing  different  amounts  of 
motion.  Figure  2  shows  mean  squared  error  vs  update  threshold.  The 
curves  of  figure  3  show  the  corresponding  compression  of  data  achieved 
with  thresholding.  The  interframe  error  curve  for  medium  motion  is 
presented  in  full  detail  in  figure  4.  Three  regions  may  be  discerned 
along  the  curve.  For  the  threshold  set  less  than  about  1.0,  all  sub¬ 
blocks  which  contain  motion  will  be  updated  together  with  some  of  the 
sub-blocks  which  do  not  contain  motion,  depending  upon  the  magnitude 
of  their  noise  variance.  For  a  threshold  set  over  the  range  1.0  to  4.0 
approximately,  all  sub-blocks  which  contain  motion  will  be  updated  except 
those  sub-blocks  which  contain  very  little  motion  (for  example,  ^ast  a 
few  pixels  in  one  corner  have  changed)  plus  those  sub-blocks  which  contain 
no  motion,  but  have  the  largest  noise  variance.  For  a  threshold  greater 
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Figure  3.2-3.  Graph  of  Interframe  Sample  Reduction 
vs.  Update  Threshold. 

1  ■  very  fast  motion,  2  =  fast  motion, 

3  =  medium  motion,  4  =  slow  motion, 

5  =  very  slow  motion. 
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Figure  3.2-4.  Graph  of  Mean  Squared  Error  vs.  Threshold, 
Medium  Motion. 
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than  about  4.0,  only  those  sub-blocks  containing  the  greatest  amount 
of  motion  will  be  updated.  This  region  of  the  wave  is  approximated 
by  the  straight  line  having  the  equation 

e2  =  0.  74  T  -  1.  63 

2 

wher^  e  is  the  mean  squared  error  and  T  is  the  update  threshold. 

This  line  does  not  fit  the  curves  for  other  degvees  o /  motion  very  well 
but  must  be  taken  as  an  overall  measure  of  the  error  of  the  system 
assuming  an  average  medium  motion. 

The  curves  for  interframe  mean  square  error  and  bit  rate 
compression  indicate  the  behavior  of  the  system  towards  the  different 
degrees  of  motion.  The  error  curves  for  medium,  fast  and  very  fast 
motion  are  reasonably  close  to  each  other  throughout  the  range  of  thres¬ 
holds  presented  while  the  corresponding  bit  rate  compression  curves 
are  nearly  linear  and  vary  from  each  other  onlv  in  slope.  In  contrast, 
the  curves  for  slow  and  very  slow  motion  differ  in  being  nonlinear  and 
not  aligned  to  the  curves  for  faster  motion.  This  effect  is  caused  by  the 
discreteness  of  the  sub-block  conditional  updating  mechanism,  which 
becomes  prominent  when  only  a  few  sub-blocks  contain  motion  and  many 
of  those  sub-blocks  overlap  stationary  areas.  The  results  for  very 
slow  motion  are  further  modified  by  noise  which  interacts  with  the  low 
difference  signal  energy,  preventing  correct  updating. 

The  data  available  for  experiments  consisted  of  five  pairs  of 
frames  representing  a  range  of  levels  of  activity  and  a  set  of  four  frames 
containing  very  active  motion.  These  frames  were  individually  coded 
by  orthogonal  transformation  and  block  quantization.  The  Hadamard 
transform  was  used  for  the  five  pairs  of  frames  with  quantization  to  3.  0 
bits  per  pixel.  The  Slant  transform  was  used  for  the  sequence  of  four 
frames  with  rates  of  3.  0  and  1.  5  bits  per  pixel.  The  result  of  interframe 
coding  the  five  pairs  of  intraframe  coded  frames,  with  a  threshold  of  4.0 
is  shown  in  figure  5.  No  discrete  errors  are  visible  as  a  result  of  the 
interframe  coding  although  the  images  are  a  little  blurred  by  the  intra¬ 
frame  coding.  Raising  the  threshold  to  9.0,  figure  >  does  introduce 
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Figure  3. 2-6. 


Interframe  Coding  with  Threshold  of  9.0. 
Reconstruction  of  Second  Frames  from  Five  Pairs  of 
Frames  Representing  a  Range  of  Activities. 


visible  error  (although  this  is  partially  masked  by  the  darkhair  of  the 
subject).  The  discrete  error  appears  at  the  edges  of  sub-blocks  which 
were  incorrectly  updated.  Two  examples  ;n  frame  number  9123  are 
the  vertical  edge  visible  or  the  subjects  left  hair  edge  and  the  L  shaped 
rick  in  the  edge  of  the  hair  to  the  right  of  and  above  the  subject's  left 
eye.  Figure  7  shows  the  results  obtained  when  coding  the  sequence  of 
four  frames  and  compares  the  effect  of  intraframe  coding  at  3.  0  bits 
per  pixel  (part  (i))  and  1.  5  bits  per  pixel  (part  (ii)).  No  errors  are 
visible  in  the  results  of  part  (i)  but  u  small  number  may  be  discerned 
in  the  second  set  of  results.  The  reason  for  these  errors  is  that  the 
higher  intraframe  compression  to  1.  5  bits  per  pixel  lowers  the  differ¬ 
ence  signal  energy.  This  is  reflected  in  the  higher  interframe  com¬ 
pression  achieved  with  these  results.  A  slightly  lower  threshold  would 
alleviate  the  problem.  The  results  also  show  that  for  the  limited 
sequence  of  four  frames  no  noise  build  up  is  evident. 

The  limited  results  indicate  that  the  scheme  is  an  effective  way 
of  coding  moving  pictures  -  the  transmission  rate  is  greatly  reduced 
while  the  coder  offers  economy  over  conventional  interframe  coders.  It 
is  not  possible  to  predict  the  effect  of  the  sub-block  update  mechanism 
in  real  time  -  it  is  possible  that  the  sub-block  structure  may  become 
visible.  Lack  of  suitable  data  has  precluded  study  of  this  problem.  The 
coder  does  produce  a  non-uniform  data  rate.  Studies  of  the  output  data 
indicate  that  the  bit  stream  may  be  smoothed  by  a  buffer  or  handled  by 
buffer  sharing  techniques  with  at  least  the  same  efficiency  as  that 
obtained  with  conventional  conditional  update  coders. 

3.  3  Quantization  Error  Reduction  for  Image  Coding 
Michael  N.  Huhns 

Quantization  is  the  process  of  representing  continuously  varying 
quantities  by  '.iscrete  intervals.  This  process  is  nonlinear  and  some  of 
the  information  about  the  original  data  is  irretrievably  lost.  The  usual 
restoration  procedure  is  to  choose  the  midpoints  of  each  quantization 
interval  as  the  estimated  values  of  the  original  data.  However  if  it  is 
known  that  the  original  data  are  correlated  and  are  non-uniformly  dis- 
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i)  Intrafj  ame  coding  at  3.  0  bits  per  pixel  ii)  Intraframe  coding  at  1.  5  bits  per  pixel 
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Figure  3.2-7.  Interframe  Coding  with  Threshold  of  4.0 
Four  Frame  Sequence 


tributed,  then  improved  restorations  are  possible  using  this  information. 
As  shown  in  a  previous  report  [l]  ,  minimum  mean  square  error  estimates 
of  correlated  data  require  the  solution  of  the  following  equation 


J"  2£P  (*) 

x=E{x|xeD]=  y -  (D 

/  pf x)  dx 
D  -  ~ 

where  x  is  the  n-dimensional  variable  to  be  quantized,  D  is  the  particular 
region  of  n- space  into  which  x  is  quantized,  and  p(x)  is  the  probability 
density  function  of  x.  A  partial  solution  to  this  equation  hi.  been  obtained 
for  data  which  have  a  jointly  gaussian  probability  distribution.  This 
solution  has  now  been  applied  to  the  restoration  of  quantized  one-dimensional 
random  signals  and  two-dimensional  transform  domain  zonal  quantized 
images.  The  results  reveal  a  decrease  in  mean  square  error  in  all  cases. 
However,  in  spite  of  the  error  reduction,  some  images  exhibit  a  degrad¬ 
ation  in  subjective  quality  after  restoration.  Hence  a  nonlinear  error 
criterion  based  on  the  human  visual  system  and  derived  by  Mannos  and 
Sakrison  [2]has  been  used  in  place  of  the  mean  square  error  function. 

Under  this  criterion  a  subjective  image  improvement  as  well  as  a 
numerical  error  reduction  are  obtained. 

To  demonstrate  the  utility  of  this  restoration  procedure,  a  randomly 
generated  gaussian  Markov  signal  has  been  quantized  and  restored.  The 
results  are  shewn  in  figure  1.  A  two  bit  per  sample  Max  quantization 
scheme  is  employed  to  obtain  the  quantized  approximation  to  the  original 
signal.  Using  this  quantized  signal  and  the  statistical  knowledge  about 
the  original  signal  as  inputs  to  the  nonlinear  estimator,  the  restoration 
decreases  the  mean  square  error  by  33%.  The  average  improvement  in 
mean  square  error  as  a  function  of  quantizing  bit  assignment  for  different 
correlation  coefficients  is  shown  in  figure  2.  It  can  be  seen  from  this 
graph  that, as  the  amount  of  correlation  in  the  Markov  process  approaches 
zero,  then  the  restoration  provides  no  error  improvement.  There  is  also 
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Figure  3.  3-2.  M.S.  E.  Improvement  for  Ouantized  Markov  Signal. 


no  improvement  as  the  number  of  quantizing  bits  becomes  large  and  the 
differences  between  the  original  signal  and  the  quantized  signal  vanish. 
Thus  the  restoration  procedure  represents  a  viable  restoration  technique 
when  the  number  of  quantizing  bits  are  small  and  the  input  signals  to  the 
quantiser  are  correlated. 

These  conditions  are  satisfied  by  the  zonal  transform  coding 
technique  for  images.  In  this  method  the  transform  samples  have  a 
gaussian  distribution:  each  is  the  sum  of  a  large  number  of  random 
variables  so  that  the  central  limit  theorem  can  be  invoked.  Now  these 
transform  samples  are  typically  quantized  according  to  a  bit  assignment 
such  as  the  one  shown  below. 


8655444433333333 

5433222222222222 

4322111111111111 
4322111111111111 
42111 11100000000 
4211111100000000 
4211111100000000 
4211111100000000 
3211000000000000 
321 1000000000000 
3211000000000000 
321  1000000000000 
3211000000000000 
3211000000000000 
321 1000000000000 
321 1000000000000 


The  restoration  technique  must  be  applied  recursively  to  these  samples 
since  it  is  only  capable  of  restoring  one  sample  at  a  time.  The  current 
best  estimates  of  the  remaining  samples  are  then  used  to  obtain  the 
estimate  of  the  sample  being  restored.  This  is  repeated  for  each  of  the 
samples  in  turn  and  then  the  entire  procedure  is  repeated  until  the 
estimates  converge.  This  convergence  has  been  found  to  be  rapid  and, 
in  most  cases,  one  iteration  is  sufficient. 

The  above  procedure  has  been  applied  to  the  images  in  figure  3  b 
and  3d.  Figure  3b  has  been  coded  with  an  average  of  one  bit  per  pixel 


by  using  a  Haar  transform  in  16  x  16  blocks,  a  zonal  coding  bit 
assignment  and  a  Max  quantizer.  For  the  quantization  and  subsequent 
restoration,  the  original  image  samples  are  assumed  to  arise  from  a 
Markov  source.  Figure  3c  is  the  restored  version  of  this  image,  uMliz- 
mg  one  iteration  of  the  estimation  procedure.  The  mean  square  error 
is  reduced  by  10%  as  a  result  of  the  restoration.  Figures  3d  and  3e, 
respectively,  have  been  quantized  to  0.  5  bits  and  restored  by  means  of  the 
above  technique.  A  reduction  of  19%  in  mean  square  error  is  obtained  in 
this  case.  Subjectively,  th*  restored  images  appear  much  less  noisy 
than  the  quantized  images  but  more  blurred,  as  is  very  evident  in  com¬ 
paring  figures  4b  and  4c.  Hence  an  error  measure  is  required  in  which 
numerical  results  match  subjective  results. 

This  has  been  provided  by  modeling  the  error  measure  after  the 
human  visual  system.  It  has  been  found  that  the  human  visual  system  is 
sensitive  to  approximately  the  cube  root  of  incident  light  intensities.  It 
is  also  most  sensitive  to  middle  spatial  frequencies  near  eight  cycles  per 
degree.  Hence  to  apply  this  error  measure,  an  image  is  processed 
according  to  the  block  diagram  in  figure  5.  The  (i,  j)th  component  of  the 
filter  function  is  chosen  to  be 

T.j  =  (.  05  +  .  11525  r)  exp  {  -(.  07125  r)1*  1 }  (3) 

where 


Figures  4d  and  4e  show  the  results  of  this  procedure  for  a  Hadamard 
transform,  with  and  without  the  restoration  step,  respectively.  r>  nere 
are  both  subjective  and  numerical  error  improvements  after  restoration. 

The  restoration  process  has  also  been  applied  to  a  color  coding 
experiment.  In  this  experiment  a  color  image  is  transformed  :o  the 
YIQ  coordinate  system  and  then  quantized  according  to  the  bit  assignment 
indicated  on  the  next  page  for  a  typical  block  of  four  pixels. 


(a)  O  riginal  Image  8  bits/pixel 


(b)  Zonal  Haar  Quantized 


(c)  Zonal  Haar  Restored 


(d)  Visual  Hadamard  Quantized 


(e)  Visual  Hadamard  Restored 


Figure  3.  3-4.  Restoration  of  0.  5  bit/pixel  transform  coded  images  for 
two  different  error  criteria. 
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Each  pixel  is  hence  coded  with  an  average  of  nine  bits.  Quantization 
restoration  provides  a  decrease  of  42%  in  mean  square  error  in  this 
case  and  an  improvement  in  subjective  quality. 

Thus  data  which  are  correlated  and  which  have  been  coarsely 
quantized  are  amenable  to  being  restored  by  the  techniques  outlined  above. 
By  choosing  a  suitable  error  criterion,  zonal  transform  coded  images 
can  be  subjectively  and  analytically  improved  so  that  they  more  faithfully 
reproduce  the  details  of  an  original  image.  Further  work  is  expected 
to  extend  the  restoration  technique  to  DPCM  and  delta  modulation  coded 
image  s. 

Reference  s 
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Criterion  on  the  Encoding  of  Images,"  IEEE  Transactions  on 
Information  Theory,  Vol.  IT-20,  July  1974,  pp.  525-536. 
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4.  Image  Restoration  and  Enhancement 


Image  restoration  and  image  enhancement  are  two  classifications 
of  image  improvement  methods.  Image  restoration  techniques  seek  to 
reconstruct  or  recreate  an  image  to  the  form  it  would  have  had  if  it  had 
not  been  degraded  by  some  physical  imaging  system.  Image  enhance¬ 
ment  techniques  have  two  major  purposes:  improvement  in  the  visual 
quality  of  a  picture  to  a  human  viewer;  and  manipulation  of  a  picture  for 
more  efficient  processing  and  data  extraction  by  a  machine.  Research 
in  both  areas  during  the  past  six  months  is  described  below. 

The  first  report  addresses  the  general  problem  of  restoration  for 
space  variant  aberrations.  Methods  are  presented  for  astigmatism  res¬ 
toration  and  restoration  methods  for  other  types  of  degradation  are  out¬ 
lined. 

The  three  reports  following  are  concerned  with  pseudoinverse 
techniques  of  image  restoration.  In  the  first  of  the  series  a  method  is 
presented  for  performing  restoration  in  the  eigenspace  of  the  degradation 
in  order  to  avoid  computational  errors.  The  second  paper  describes 
a  restoration  technique  which  seeks  to  ensure  that  the  estimated  image 
pixel  values  are  neither  negative  nor  greater  than  an  upper  bound.  In 
the  third  paper  a  computational  algorithm  for  pseudoinveroe  restoration 
is  developed.  With  this  algorithm  computation  is  reduced  significantly 
compared  to  conventional  methods  and  numerical  stability  is  also  improved. 

The  next  paper  discusses  the  use  of  spline  function  interpolation 
in  the  development  of  an  image  restoration  algorithm,  and  in  the  next 
paper  the  effects  of  non-uniform  sampling  of  a  blurred  image  are  analyzed. 
It  is  shown  that  proper  non-uniform  sampling  actually  can  improve  the  qual¬ 
ity  of  an  image  restoration. 

Finally,  a  method  of  histogram  exponentiation  for  image  detail 
enhancement  is  described.  With  this  process  a  grey  level  mapping  is 
performed  to  produce  an  output  image  whose  histogram  follows  an 
exponential  density. 
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4.  1  Restoration  for  Space- Variant  Aberrations 

Alexander  A.  Sawchuk  and  M.  Javad  Peyrovian 


Several  previous  publications  have  detailed  the  difficulty  of 
achieving  image  restoration  for  the  most  general  case  when  the  linear 
degrading  system  is  space- variant.  By  systematically  examining  the 
sources  of  degradation  and  using  all  available  a  priori  knowledge,  con¬ 
siderable  simplificatic  of  the  problem  can  be  achieved  [  1,  5]  . 

An  example  of  restoration  for  images  of  100  x  100  pixels  degraded 
by  third -order  astigmatic  aberrations  has  been  previously  described, 
and  descriptions  of  space- variant  point- spread  functions  (SVPSF)  have 
been  given  [  1,  5]  .  Recent  progress  in  this  area  has  permitted  images 
with  astigmatism  of  up  to  128  x  128  pixels  in  size  to  be  restored  on  the 
USC  360/44  computer  system  using  an  SVD  technique.  With  the 
addition  of  memory  and  a  faster  processor  now  underway,  these  capa¬ 
bilities  should  be  increased  in  the  near  future.  Another  improvement  in 
the  simulation  and  testing  of  the  astigmatism  algorithm  has  resulted 
from  the  use  of  an  improved  quadrature  formula  to  evaluate  the  space- 
variant  integral 


J,xl’x2)  '/ /h,’VV'Yu 2>°<ui-u2>  duldu2 


(1) 


describing  the  degradation  process.  Here  ^(u^u^)  and  -P(Xj,x^)  are 
the  object  and  image,  respectively,  and  h(x  ,  x  ,u  ,u  )  is  the  degrading 

lulu 

SVPSF.  The  new  quadrature  formula  involves  the  expansion  of 
h(x  ,x  ,u  ,u_)  into  spline  functions,  thus  eliminating  artifacts  in  the 

i  ^  1  w 

digital  simulation  of  the  imaging  system  and  providing  a  better  test 
of  the  restoration  method. 

A  test  of  space-variant  restoration  for  the  most  difficult 
third-order  aberrations  of  combined  astigmatism  and  curvature  of  field 
is  now  in  process.  The  theory  behind  the  technique  is  an  extension  of 
the  method  for  astigmatism  alone.  The  general  degradation  of  eq.  (1) 
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is  first  rewritten  with  a  polar  coordinate  transformation  in  the  form 


J<vV  -If  h(xr,  ur,  xp.  u0)  uQ)  durduf 


(2) 


where  (x  ,  x  )  and  u  ,  u  )  are  image  and  object  polar  coordinate  vari- 
r  9  r  p 

ables.  xnis  operation  is  digitally  performed  on  the  degraded  image  by 
rearranging  samples  in  memory.  A  simplification  results  because  of 
the  circular  symmetry  of  many  optical  systems.  In  this  case,  the  system 
aberration  function  are 


x  -u  =  (2C+D)u  e  cos  e 
r  r  r  r 


(3a) 


X--U-  =tan  *[Du  e  sin  e  /l  +  (2C+D)  u  e  cos  e  ]  =g(u  )  (3b) 

OH*  rr  n  rr  o  r 

where  C  and  D  are  the  astigmatism  and  curvature  of  field  coefficients 
respectively,  and  e  and  are  exit  pupil  variables.  The  SVPSF' s 
derived  from  eq.  (3)  and  limiting  cases  of  C  =  0  and  D  =0  can  be  easily 
derived  [  l]  . 

By  observing  that  the  radial  blur  of  eq.  (3a)  is  a  function  of  radial 
coordinates  only,  and  that  the  8  blur  of  eq.  (3b)  is  a  slowly  varying 
function  g(ur),  the  next  step  in  restoration  is  to  rewrite  eq.  (2)  as 

00 

J(Vxrf  =/ /  h(xr.Ur,x0.uB,g(ur))O(ur,ufi)  durdu0  (4) 

to  emphasize  the  dependence.  If  D  =0  in  eq.  (3)  the  degradation  is  pure 
astigmatism  without  curvature  uf  field  and  the  blur  is  entirely  radial 
with  no  blurring  in  9  . 

Defining  a  Fourier  transform  of  J(x^,Xq)  in  the  x^  variable  by 

00 

J(xr>\)  =  J ^<xr.x9)  exp(-j2rrXxg  )  dx0  (5) 
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the  transform  of  both  sides  of  eq.  (4)  is  taken  to  obtain 


J  (x  ,  X)  =  /  /  ^(u  ,  u  )h(x  ,u  ,X,g(u  ))  exp(-j2nXu  )du  du 
r  j  j  rn  rr  r  bar 


_  00 


(6) 


where  h  is  the  transform  of  h  in  xfl.  Grouping  terms  containing  u 

b  8 

on  the  right  side  of  eq.  (6)  enables  a  transform  in  this  variable  to 

be  evaluated.  The  resulting  transformed  function  (3(u  ,  X)  is  given 

r 

by 


0( 


ur>  M  =  j 0(ur.ue)exp(-j2nXue)du( 


(7) 


and  the  reduced  system  equatiin  obtained  from  eq.  (4)  is 


«^(x  ,  X)  -f  h(x  ,  u  ,  X)  <3(u  ,  X)du 
r  J  r  t  r  t 

_  00 


(8) 


where 


h(x  ,u  ,  X)  =  h(x  ,  u  ,  X,  g(u  )) 
r  r  r  r  r 


(9) 


is  rewritten  as  a  function  of  three  variables  to  show  explicit  dependence. 

This  derivation  shows  the  procedure  to  be  used  for  the  astig¬ 
matism  and  curvature  of  field  restoration.  Following  a  polar  coordinate 
transformation,  a  Fourier  transform  in  x  as  expressed  by  eq.  (5)  is 
performed  to  partially  decouple  8  blur  as  a  slowly  varying  function  of 
giu^).  The  reduced  system  is  then  given  by  eq.  (8),  and  an  estimate 
of  <3(ur>X)  is  produced  by  singular-value  decomposition  (SVD)  techniques 
[  l]  for  each  separate  X  by  techniques  similar  to  that  used  for  astig¬ 
matism  alone.  Efforts  will  be  made  to  reduce  the  computational  effort 

required  in  this  part  by  using  the  known  variation  of  h(x  ,u  ,  >)  with  X. 

_  r  r 

After  the  entire  ^(u^,  X)  has  been  obtained,  a  .  eries  of  one-dimensional 

inverse  transforms  in  x  is  taken  to  find  0(u  ,u„),  and  an  inverse  polar 

r  0  r 

coordinate  distortion  is  used  to  get  <3(u  ,  u  )  as  the  final  restored  object. 

X  Lt 
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This  procedure,  while  requiring  large  capabilities  in  computing  and 
storage,  is  the  only  method  for  doing  restoration  on  images  of  even 
moderate  size.  The  general  four-dimensional  space-variant  blur  of 
eq.  (1)  is  effectively  reduced  to  a  set  of  space- variant  two-dimensional 
problems  whose  point- spread  functions  depend  in  a  well-behaved  way  on 
X  . 

Following  tests  of  this  procedure,  future  work  on  this  problem 
will  include  modifications  for  noisy  data  and  an  analysis  of  sampling 
problems  and  approximations  in  performing  a  continuous  space-variant 
restoration  by  computer.  The  analysis  as  presented  is  for  a  general 
aperture,  although  additional  simplification  may  be  possible  for  other 
specialized  aperture  shapes. 
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4.  2  Image  Restoration  in  the  Eigenspace  of  the  Degradation 


Harry  C.  Andrews  and  Monty  Adler 

It  has  recently  been  noted  that  there  are  certain  mathematical 
properties  associated  with  point  spread  function  blur  matrices  that 
make  digital  image  restoration  extremely  attractive  in  the  eigenspace 
of  the  degradation  [  l]  .  Consider  the  separable  space  variant  imaging 
model  of 


G=AFB  (1) 

where  G  is  the  matrix  of  available  image  data  points  and  Fis  the 
array  of  object  data  points  which  are  of  interest  for  estimation.  Usually 
both  A  and  B,  the  column  and  row  blur  matrices  respectively,  are 
nearly  singular.  Penrose  [  2 ]  investigated  this  model,  without  regard 
to  imaging  systems,  and  suggested  the  pseudoinverse  restoration  which 
provides  an  estimate  of  the  object  as 


~  +  + 

F  =  A  G  B  (2) 

where  _A  and  B  are  Moore- Penrose  pseudoinverses  respectively 
[  3] .  The  blur  matrices  A  and  B  are  componentwise  positive  and 
experimentally  appear  to  be  oscillatory  [  4,  5]  ,  The  oscillatory  property 
implies  that  the  eigenvectors  associated  with  decreasingly  smaller 
eignvalues  have  increasingly  larger  zero  crossings.  Examination  of  the 
pseudoinverse  solution  of  eq.  (2)  in  eigenspace,  reveals  that  the  solution 
can  be  cast  in  the  form 


~  +  T  T  T 

F  =  V  A  U  GV  A  U, 
a  "a  a  b— b  ~  b 


(3) 


where 


A  =  U  A  V 
—  ~  a~a—  a 


(4a) 


and 


B=U.^E 


T 

Lrlr^b 


(4b) 


All  matrices  are  orthogonal  with  A  and  _A^  being  diagonal.  The 

+  -f 

matrices  A  and  A,  are  also  diagonal  with  non-zero  entries  being  re- 
— a  — b 

ciprocals  of  the  entries  of _A^  andj^  resPectively*  except  for  zeros  placed 
on  the  diagonals  of _A*  and  J\^  for  small  (approximately  zero)  diagonal 
entries  of _A&  and_A^.  By  noting  that 


~  T  T  T 
F  =  V  A  a  A  U, 

—  —  a~a — b  — b 


(5) 


where 

a  =  U  TGV  (6) 

-  -a  — b 

it  is  found  that  the  a  matrix  is  the  image  G  represented  in  the  eigenspace 
of  A  and  B.  Because  of  the  oscillatory  properties  of  A  and  B,  con¬ 
tributions  to  F  associated  with  large  X  1  coefficients  have  correspond¬ 
ingly  large  zero  crossing  eigenvectors.  These  eigenvectors  are  then 
correlated  with  the  image  G  which  tends  to  be  slowly  varying;  and  as 
such,  the  correlation  tends  to  zero.  But  those  eigenvalues  of  least  con¬ 
fidence,  due  to  uncorrelated  noise  and  computational  error,  are 
associated  with  eigenimages  of  many  sign  changes,  which  in  turn  correlate 
to  zero  with  G.  Consequently  one  encounters  the  happy  circumstance 
of  computational  error  not  being  as  critical  as  could  be  as  a  result  of 
effective  zero  correlation  of  the  image  G  with  the  corresponding  com¬ 
putationally  questionable  eigenvalued  **igenimages.  Thus  greater  success 
can  be  anticipated  with  this  pseudoinverse  approach  for  image  restoration 
than  for  restoration  of  signals  with  correspondingly  higher  frequency 
content.  Restoration  algorithms  based  upon  these  observations  are 
currently  under  investigation. 
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4.3  Pseudoinverse  Method  of  Bounded  Image  Restoration 
Harry  C.  Andrews  and  Monty  Adler 

The  space  variant  separable  imaging  equation  given  by 

G  ---  A  F  B  (1) 

has  been  investigated  from  a  restoration  viewpoint  in  previous  reports. 
With  A  and  B  being  the  column  and  row  blurs  of  the  object  F  respectively, 
the  image  G  will  be  used  as  a  means  of  estimating  F.  Traditional 
solutions  to  this  problem  utilize  the  pseudoinverse  in  which 

A  +  + 

F  =A  GB  (2) 

,  .  +  + 

and  A  and  B  are  pseudoinverses  of  A  and  B  respectively.  Usually  the 
blur  matrices  are  nearly  singular.  Unfortunately,  although  the 
pseudoinverse  solution  provides  the  best  estimate  of  F  in  a  least  squares 

A  ^  ^ 

sense  for  a  minimum  norm  F,  ( i.  e.  F  minimizes  l|  F  ||  ),  the  estimate 
does  not  take  advantage  of  all  a  priori  knowledge.  For  example,  the 
solution  ignores  the  fact  that  the  elements  of  F  are  intensities,  which 
cannot  go  negative  (therefore  positive  restoration),  and  which  cannot 
exceed  some  physically  realizable  total  light  energy  (therefore  bounded 
restoration).  To  utilize  the  positive  bounded  restoration  model  the 
pseudoinverse  solution  of  eq.  (2)  is  suggested  as  an  initial  condition  to 
a  nonlinear  programming  algorithm  to  guarantee  positivity  and  upper 
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boundness. 


Towards  this  end,  consider  the  pseudoinverse  solution 


where 


and 


+  T  +  T 

F  =  V  A  U  G  V  A  U 
—  — a~a  — a - b~ b  —  b 


(3) 


A  =U  A  V 
—  — a- a  — a 


(4a) 


§  *J?b4  4 


(4b) 


The  objective  function  for  normalization  purposes  will  be  to  minimize 
G-G  ^  subject  to  the  bounded  restoration  constraint  on  F.  Thus 


W  =  II  G-G  I!  2 


(5a) 


W  -  !!  A(F -F)  B  ||  (5b) 

w  *  ||a||2  ||f-f||2  ||bI|2  (6) 

Therefore,  for  analysis  purposes  one  can  choose  to  minimize  ||  F-F  l|  . 
However,  since  the  pseudoinverse  is  simply  being  used  as  an  initial 
condition  and  since  the  eigenspace  of  the  degradation  is  a  useful  domain 
for  restoration  processing,  the  variables  of  optimization  are  restricted 
to  be  diagonal  in  the  eigenspace  domain.  Thus  the  estimate 

*  +  T  +  T 

F  =  V  6  A  U  G  VJ,  &  U.  (7) 

—  —  a  ~a  ~a  ~a - b— b  — b  —  b 

is  found  by  adjusting  the  ^  diagonal  matrices  such  that  ||F-F  l|  ^  is 

a 

minimized  and  £  is  componentwise  positive  and  bounded.  This  reduces 

2 

the  restoration  problem  to  a  2K  variable  problem  with  N  positive 
bounded  boundary  constraints.  Here  K  is  the  number  of  non  zero 
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eigenvalues  retained  in  the  F  estimate.  Since  the  solution  is  in  the 
eigenspace  domain,  the  effect  of  each  6^  variable  is  felt  throughout  the 
entire  estimated  object  F,  and  as  such  is  much  more  effective  than 

A 

simply  constraining  a  component  entry  (pixel)  of  F  itself.  A  Fiacco 
and  McCormick  [  l]  programming  algorithm  has  been  utilized  to  adjust 
the  weights  in  degradation  eigenspace  to  obtain  positive  bounded 
restorations  with  some  success.  While  this  work  is  still  in  progress, 
it  appears  that  convergence  is  very  rapid  due  to  the  pseudoinverse 
starting  point  whereby  it  is  meant  that  an  iterative  pseudoinversion 
technique  is  developed  to  obtain  F  until  the  positive  bounded  con- 
straints  are  violated  at  which  time  the  nonlinear  programming  algorithm 
is  call  to  correct  for  the  near  singularity  in  the  K**1  pseudoinverse 
solution. 
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4.4  A  Fast  Pseudoinverse  Image  Restoration  Algorithm 
William  K.  Pratt  and  Faramarz  Davarian 

It  is  often  possible  to  model  an  image  degradation  process  by  the 
vector  equation 

J  =  B  £+  n_  (1) 

where  g  denotes  a  column  scanned  M  x  1  vector  of  physical  samples  of 
the  blurred  image;  Ms  an  N  x  1  vector  of  column  scanned  points  in  the 
ideal  image  field,  B  is  an  M  x  N  blur  operator  matrix  representing 
a  convolutional  blur;  and  ii  denotes  an  M  x  1  vector  of  observation  noise 
or  uncertainty.  For  this  model,  the  ideal  image  vector  can  be  estimated 
by  pre-multiplication  of  the  observation  by  the  generalized  inverse, 

B  ,  of  the  blur  matrix.  Thus, 
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(2) 


f  =  B  g 


where  B  is  an  N  x  M  matrix  that  can  be  computed  by 


T  -1  T 
B  =  (B  B)  B 


(3a) 


if  B  is  of  full  column  rank,  or  can  be  computed  by 


B 


T  -1 
(B  B  ) 


(3b) 


if  B  is  of  full  row  rank.  The  solution  of  eq.  (2)  is  a  minimum  mean 
square  error,  minimum  norm  estimate.  There  are  two  major  difficulties 
with  pseudoinverse  restoration:  if  the  noise  level  is  high,  the  solution 
may  be  unstable  as  a  result  of  the  usual  ill  conditioning  of  the  blur 
matrix:  and,  computation  of  the  generalized  inverse  and  restoration  by 
eq.  (2)  is  usually  a  large  task.  The  former  problem  can  be  avoided  by 
restoration  constraints;  consideration  is  given  here  to  efficient 
computational  techniques  for  pseudoinverse  restoration. 

Fast  Pseudoinverse  Algorithm.  As  a  simplification  in  the  develop 
ment  of  the  algorithm,  consideration  will  be  initially  limited  to  a  one 
dimension'll  model  in  which  the  ideal  image  is  represented  by  quadrature 
points  at  its  Wyquist  rate  and  the  blurred  image  is  sampled  at  the  same 
rate.  Then,  let  the  impulse  response  be  represented  by  the  Lxl 
vector  h.  The  blur  matrix  of  eq.  (1)  then  assumes  the  form 


h(L)  .  .  .  h(l)  0  .  .  .  0 
0  h(L)  .  .  .  h ( 1)  ...  0 


B  = 


0  h(L)  .  .  .  h(  1) 


(4) 


In  this  case  the  number  of  ideal  image  points  and  the  number  of  observed 
image  samples  are  related  by 
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N  =  M  +  L-l 


(5) 


Now,  let  two  vectors  and  g^,  be  formed  by  selecting  the  center 
portions  of_f  and  g,  respectively.  These  truncated  vectors  are  obtained 
by  dropping  L-l  elements  at  each  end  of  the  appropriate  vector  by  the 
operations 


=S2(K) 
— N 


K  =  N-  2(L- 1) 


It  =  — 


(R) 

M 


R  =M-2(L-1) 


where 


L-l  K  L-l 


for  J  =  K  +  2(L-1).  Figure  la  illustrates  the  relationship  of  all  vectors 
for  N  =  9  original  vector  points,  M  =  7  observations,  and  an  impulse 
response  of  length  L  =  3. 

Suppose,  now  that  the  data  sequence  f^  is  discretely  convolved 
with  the  impulse  response  sequence  n  yielding  the  output  sequence  q 
as  defined  by  the  vector  equation 

i  =  P  iT  (6) 
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(a)  sampled  continuous  convolution 


(b)  discrete  convolution 


Figure  4.4-1.  Examples  of  one  dimensional  sampled  continuous 
convolution  and  discrete  convolution. 


where 


Figure  lb  illustrates  the  relationships  between  the  vector  sequence  for 
discrete  convolution.  An  estimate  of  f^  can  be  obtained  by  premultiplying 
the  output  vector  c^by  the  generalized  inverse  of  D.  That  is 

iT  =  D'a  (8) 

Referring  to  figure  1,  it  is  observed  that  the  elements  of  g^,  are 
identical  to  the  center  elements  of  q.  Thus, 


tsE’Ji-isSN 


while  the  remaining  end  elements  differ  in  general.  Now,  let  a  matrix 
W  be  defined  which  operates  on  the  physical  sample  vector  g  to  produce 
an  approximation  to  the  output  vector  q  for  discrete  convolution. 


The  structure  of  W  will  be  derived  later.  Then,  an  estimate  of_f,j, 
is  formed  by 
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no) 


It 


D  q  =D  W  g 


By  this  procedure  an  estimate  of  the  center  part  of  the  image  vector 
f  can  be  obtained  by  use  of  the  generalized  inverse  operator  D  rather 
than  the  generalized  inverse  operator  B  .  The  advantages  of  the 
former  procedure  over  the  latter  are  that,  in  the  absence  of  noise,  the 
solution  obtained  with  D  is  unique  and  exact,  whereas  there  are  an 
infinite  number  of  feasible  solutions  for  the  B  model.  Also,  with  the 
D  operator,  it  is  possible  to  perform  the  restortion  by  Fourier  domain 
processing  quite  efficiently. 

Consideration  will  now  be  given  to  the  structure  of  the  weighting 
matrix  W.  The  objective  of  weighting  is  to  express  the  vector  q  in  terms 
of  the  elements  of  the  observation  g.  For  the  example  of  figure  1, 

q(l)  =  g(l)  -  h( 3)f ( 1)  -  h(2)f(2) 
q ( 2)  =  g(2)  -  h(3)f(2) 

q(m)  -  ^(m)  2<m<M-2  (11) 

q(M-l)  =  g(M-l)  -  h(l)f(N-l) 

q(M)  =  g(M)  -  h(2)f(N-l)  -  h(l)  f(N) 

Since  the  values  of  £are  not  known,  the  correspondence  of  eq.  (11) 
cannot  be  made  directly.  However,  by  making  an  assumption  on  the 
continuity  of  the  original  image  vector  that 


and 


f(l)  =  f  ( 2)  =  f(3) 


f(N-  2)  =  f(N-l)  =  f(N) 


then  it  is  found  that 
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5,1, .  Mi 

q(2,  =  g(2)[h(D  +  h(2>] 

s 


where  S=  h(l)  4-  h(2)  +  h(3)#  Similar  equations  exist  for  q(M-l)  and 
q(M).  This  procedure  can  be  generalized  for  any  size  vectors.  Also, 
more  complex  prediction  algorithms  may  be  employed.  Figure  2 
illustrates  the  expected  mean  square  restoration  error  ofj:^  for 
various  prediction  algorithms  as  a  function  of  the  correlation  of  elements 
of  f_  under  the  assumption  that  {_  is  a  sample  of  a  Markov  process  with 
correlation  factor  p.  A  first  order  predictor  provides  a  significant 
improvement  over  the  zero  order  predictor,  and  for  high  correlation 
factors  even  surpasses  the  estimate  for  the  underdetermined  model. 

The  efficiency  of  the  computational  algorithm  is  based  upon  the 
use  of  a  circulant  blur  operator  defined  as 


h(l)  0 
h(2)  h(l) 
h(2) 

•  • 

h(L)  h(L-l) 

ML) 

•  • 

0  o 

0  0 


0  h(3)  h(2) 

•  • 

•  • 

h(L)  . 

0  h(L) 


0  h(l)  0 

0  h(2)  h(l) 


Then,  for  the  noise  free  case,  it  can  be  shown  that 
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And  hence,  the  image  estimate  may  be  obtained  by 


R 

i 


J-R 


where  C_  is  the  generalized  inver se  of  C,  which  i s  al so  circulant.  Since, 

the  two  dimensional  Fourier  transform  of  a  circulant  matrix  is  of 

diagonal  form,  it  becomes  more  efficient  to  perform  the  computations 

2 

in  the  Fourier  domain.  The  J  matrix  operations  associated  with  the 
pseudoinverse  multiplication  can  then  be  replaced  by  J  scalar  multiplies 
plus  2  J  log2J  operations  required  for  the  fast  Fourier  transformations. 
While  the  fast  pseudoinverse  algorithm  generalizes  quite  easily  for  two 
dimensional  image  fields.  Application  of  the  technique  to  image  res¬ 
toration  is  now  underway. 


4.  5  Spline  Function  Image  Restoration. 

Steve  Hou  and  Harry  C.  Andrews 

For  an  estimated  object  the  continuous  -  discrete  imaging 
system  model  can  be  represented  as  the  discrete  matrix  equation 

(A  +  yB)  £  =  d  (1) 


or  as 


A  _  A  A 

(ATA  +  YB  X  B2)  c  =  ATg 


using  cubic  B  spline  interpolation  for  the  estimated  object, 
eqs.  (1)  and  (2) 
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The  size  of  the  image  is  I  x  J  and  that  of  the  interpolating  grids  for  the 
estimate  object  is  M  x  N.  The  term  S  (C)  denotes  the  one-dimensional 

•*  K 

cubic  B  spline  function  centered  at  grid  point  x^  and  d^  and  £  are 

column  vectors  of  the  image  g. .  with  unknown  coefficients  c  arranged 

ij  Kv 

in  a  lexicographic  order.  The  symbol  ®designates  a  Kronecker  (tensor 
or  direct)  product  between  two  matrices  [  l]  .  Both  matrices 
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and 


T 

A  =  A  A 

B  =  B,  ®  B 
-  -1  -2 

are  real,  symmetric,  non-negative  definite,  and  the  diagonal  terms 
are  larger  than  any  off-diagonal  terms.  Furthermore  the  entries  of  A 
are  positive.  Making  use  of  the  following  equal  energy  constraint 
on  the  imaging  system 

I  J 

E  E  vc>M- 1 

i=l  j=l 

for  all  C  and  H  ,  A  can  be  shown  to  be  a  Markov  matrix  [  l]  . 

Ultilizing  the  special  properties  of  cubic  B  spline  function  one 
can  further  show  that 


Thus  Bj  is  a  banded  matrix  of  seven  element  wide.  Also,  B{  is  a 

Toeplitz,  almost  cyclic  and  strictly  diagonally  dominant  matrix  with 

the  properties  mentioned  previously.  The  matrix  f^has  the  same 

entries  as  B  except  they  may  have  different  size.  Hence  B  and  B 
1  “1-2 
are  positive  definite  matrices  C  2]  . 

Iterative  methods  have  been  formulated  for  solving  the  unknown 

coefficients  c^  in  eq.  (2).  The  matrices  A  and  B  usually  have  very  large 

dimensions  in  an  image  restoration  problem,  hence  one  is  fa  ed  with 
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the  problem  of  inverting  a  matrix  of  very  large  dimension.  It  is 
known  that  the  number  of  multiplications  in  inverting  matrices  grows  as 
a  nonlinear  function  of  the  dimensions  of  the  matrices  being  inverted. 
Therefore,  for  large  dimension,  matrix  inversion  is  a  non-trivial  operation 
even  with  a  modern  digital  computer.  For  this  reason  it  is  assumed 
that  the  point  spread  function  is  separable  i.  e. 


h  .  (C.k)  =  h.(C)  h.(H)  (3) 

U  1  J 

for  all  C>Handi,j.  Now  matrix  A  becomes 


A  =  P  ®  P 

_  — x  — y 


where 
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—  X 
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-y 
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(4) 


Then  the  following  iterative  equation  can  be  derived  from  eq.  (2) 
by  making  use  of  eq.  (4). 


:(l)  =  c(0)  -  y(Aj  +  Bj)  c(i"1)(b2a2'1 


(5) 


Where  C  is  the  .....crix  of  the  unknown  coefficients  >  +  designates 
the  pseudoinverse  [3],  and 
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Work  is  in  progress  to  determine  the  unknown  parameter  y 
(o<  Y<  1 )  from  the  iteration  procedure  of  eq.  (5)  and  the  constrainted 
iterative  solution  of  C  imposed  by  the  positive  restoration  and  equal 
energy  constraints. 
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4.6  Nonuniform  Sampling  of  Observation  Space 
Faramarz  Davarian 

The  fundamental  model  describing  the  image  restoration  pro¬ 
blem  under  the  assumption  of  a  space  invariant  blur  function  is  given 
by 

OC 

fo(x)  =  f  fj(a)g(x  -  a)da  (1) 

-00 

Here,  fj(a)  represents  an  ideal  image  line  and  rjx)  an  observed  image 

line.  To  restore  fjta)  numerieally,  the  above  integral  equation  may  be 
discretized  as  follows 

N 

W  Z  CU  -  J>  (2) 

j=l 

where  f^j)  are  uniformly  spaced  nodes  of  the  quadrature  formula, 
fQ(x.)  are  samples  of  the  observation  and  c  are  quadrature  coefficients. 
Equation  (2)  can  be  represented  in  vector  form  as 
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(2) 


where 

B. .  =  c. .g(x.  -  j  ) 

»J  U  ‘ 

The  problem  of  image  restoration  has  now  been  reduced  to  a 
regression  problem;  given  the  observed  vector _f  and  the  blur  matrix 
B,  a  suitable  e.Hmate  °f_fj  must  be  found. 

If  in  eq.  (1)  matrix  B  has  full  column  rank  (M^  N),  the  model 
is  called  overdetermined.  Under  the  assumption  of  full  column  rank 
the  pseudoinverse  of  B  is  defined  to  be  [  l] 


+  T  -1 
B  =  (B  B) 


and  the  estimate  of_f^  is  given  by 


Condition  Number.  The  condition  number  of  a  system  is  a  measure 
of  the  affect  of  an  input  perturbation  (input  noise)  on  the  output  of  the 
system.  Assume  that  the  observed  image  function  has  been  perturbed 
by  observation  noise  by  an  amount  of  A_f^.  The  error  in  the  estimate  of 
f  is  bounded  by  [  2  ] 


(3) 


where  C 

a 


is  the  condition  number  of  B  and  is  given  by 


C 


B 


=  I|b||  . 
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From  inequality  (3)  the  importance  of  the  condition  number  becomes  clear. 
A  large  condition  number  results  in  numerical  unstability;  a  small 
observation  error  will  cause  a  large  error  in  the  estimation  of  the  ideal 
image.  A  linear  system  with  relatively  large  condition  number  is  often 
called  ill  conditioned. 

Equally  Distanced  Samples  of  Observation.  Assume  that  the 
original  image  is  sampled  at  points  1,2, . . . ,  N,  and  the  sampling 
interval  has  unit  length.  Suppose  the  sampling  period  for  the  observed 
image  is  Ax(Ax  s  1)  and  the  length  of  the  degrading  function  is  L;  g(t)  =  0, 
if  |  1 1  >  L/2.  The  relationship  between  the  number  of  observed  samples 
M  and  the  number  of  original  samples  N  is  then  given  by 


M 


N  -  L 
Ax 


+ 


1 


For  a  given  N  and  L,  M  should  be  chosen  to  minimize  the 

condition  number  for  B.  Consider  the  set  of  blur  matrices  B  for 

— M 

which  M  is  sequentially  set  equal  to  N,  N  +  1,  N  +  2,  ...  ,  etc.  and  the 
correspond'ng  condition  number  is  C__  _  .  Figure  1  contains  a  typical 
plot  of  the  condition  number  as  a  function  of  M.  Usually  the  shape  of 
such  a  condition  number  curve  depends  on  the  variance  of  the  blur 
function.  In  general,  the  curve  does  not  decrease  monotonically.  The 
curve  assumes  its  peak  when  M  =  N,  and  as  M  grows  the  curve  tends  to 
decrease  with  some  periodic  upward  jumps.  The  period  is  related  to 
Ax  and  the  fact  that  for  some  integer  n,  nA  x  =k  (in  some  cases  such  n 
may  not  exist)  where  k  is  an  integer.  The  period  of  the  jumps  is  usually 
a  function  of  n  and  N.  For  very  large  M,  the  curve  tends  to  increase,  but 
the  amplitude  of  the  jumps  decrease. 

One  method  to  improve  upon  the  condition  number  of  B  is  non¬ 
equal  spacing  of  the  observed  samples.  It  is  clear  that  the  middle  points 
of  the  ideal  image  line  contribute  to  more  points  of  observation  than  the 
points  which  are  closer  to  the  line  boundary.  If  the  space  limited 
degradation  function  exists  only  on  an  interval  of  length  L,  a  typical 
pixel  in  the  middle  would  contribute  to  all  the  points  in  the  observation 


N  2N 

NUME 


Figure  4.6-1.  Typical  shape  < 
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image  line  which  lie  in  an  interval  of  size  L.  On  the  other  hand  the 
two  samples  at  the  very  ends  of_fj  only  contribute  to  two  samples  at  the 
end  points  of_fQ.  This  simple  argument  suggests  that  non-uniform 
sampling  of  observed  image  would  bring  out  all  the  information  needed 
to  easily  reconstruct  the  original  image. 

The  ill- conditioning  associated  with  uniform  sampling  of  a  blurred 
image  have  been  examined  experimentally.  In  the  experimental  model 
a  Gaussian  blur  function  has  been  employed  with  16  equally  spaced  ideal 
image  samples  on  an  interval  of  length  15,.  The  variance  of  the 
Gaussian  blur  function  is  unity  and  its  space  limited  length  is  6  pixels. 

For  16  equally  spaced  samples  of  the  observation,  the  condition 
number  is  about  300,  000.  When  the  number  of  observation  samples  in¬ 
creases,  the  condition  number  tends  to  decrease.  For  values  of  M 
above  50  the  condition  number  stays  under  2,000  ,  but  it  never  becomes 
much  less.  As  mentioned  before  C0  is  not  a  monotone  function,  i.  e.  , 
for  some  values  of  M  it  may  increase  by  increasing  M.  For  a  suitable 
non-uniform  sampling  of  the  observation,  the  condition  number  becomes 
about  2,000.  The  different  methods  to  find  the  locations  of  the  non- 
uniform  samples  are  discussed  in  the  next  section.  Note  that  there  is 
an  advantage  in  non-uniform  sampling  with  M  =  N  as  compared  to  uniform 

sampling  with  M  s  N,  since  for  the  first  case  it  is  necessary  only  to 

.  T 

invert  matrix  B  itself  whereas  for  the  second  case,  B  B  must  be 

inverted. 

Spacing  of  Non-uniform  Samples.  The  two  end  points  of  the 
observation  line  play  an  important  role  in  the  sampling  process;  the  be¬ 
ginning  and  the  end  points  of  the  observed  line  are  essential  to  enable 
the  system  to  reconstruct  the  object.  To  simplify  the  explanation, 
consider  the  left  half  of  the  observation  only.  Starting  from  the  beginning 
of  the  left  half  line,  call  the  unit  distance  points  of  the  observation 
xl*  X2’  *  *  *  e*c»  as  shown  in  figure  2a.  In  the  same  manner,  the  points  of 
the  object  line  are  denoted  as  y  fy  ,  ...  etc.  The  problem  then  is  to 
estimate  points  y^.y^....  etc;  using  N  samples  of  the  observed  intervals. 
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observation  line 


(b)  nonuniform  sampling 


Figure  4,6-2.  Relationship  between  object  and  observation 
line  samples. 


Interval  [  x  ,  x  )  is  the  only  one  containing  information  about 

1  4 

[v  .v  )  .  The  interval  [x. ,  x  ) ,  also,  contains  information  about  the 
71  2  12 

larger  interval  [  y  ,  y  ) .  Now,  consider  the  interval  [  y^,  y^)  .  There 

are  only  two  intervals  of  the  observed  line  which  can  contribute  to 

the  restoration  of  [  y2,y3)  :  ^  x^,  x^)  an<*  ^  X2’X3^  '  In  t^ie  intervai 

[y,,yl  ,  only  [  x  ,  x  ),  [x_,  x  J  and[x  ,  x  )  of  the  observation  contain 
3  4  12  23  34 


information  about  [  y^  y  )  .  This  simple  argument  can  be  generalized 


to  a  typical  unit  size  interval  [  y^y^)*  There  are  at  most  L  -  1  of  the 


(L  is  the  length  of  the  degrading  function)  observation  intervals  which 


could  contribute  to  reconstruction  of  [  y.fy.,,).  The  number  of  intervals 

3  J+l 


relating  to  [y^.y^)  equals  j  if  j  <L  ,  and  equals  L-lifj^L-1 


(assume  y.  is  on  the  left  half  line).  Suppose  for  every  interval  of  the 
ideal  line  it  is  necessary  to  have  k  samples  of  the  observed  line.  For  the 


first  interval  [y^.y^,  all  these  k  samples  must  lie  on  [  x^,  x^).  For  the 
second  interval  [  y^.y^),  half  must  lie  on  [  x^.x^),  and  half  on  [  x^  x^). 
Likewise,  one  can  consider  the  sample  distribution  for  the  <ther 


intervals. 


Figure  2b  demonstrates  the  method.  In  the  table,  the  number 
of  samples  which  are  needed  to  reconstruct  the  corresponding  ideal 
line  interval  is  listed  before  each  observation  interval.  Here,  16 
samples  of  the  object  line  are  uniformly  placed  on  an  interval  of  size 
15,  and  the  size  of  the  degrading  function  is  6.  The  length  of  the  obser¬ 
vation  line  is,  therefore,  15  -  6  =  9.  Three  samples  are  placed  on  the 
first  interval,  two  on  the  second,  and  one  is  placed  on  the  third.  The 


same  distribution  is  considered  for  the  last  three  intervals.  The  middle 


of  the  line  is  sampled  exactly  as  the  object  line. 

Figure  3  contains  a  blurred  picture  and  its  restoration  using  the 
technique  described  in  this  report.  The  restoration  is  exact  since  it  lias 
been  performed  in  the  absence  of  noise. 
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4.7  Histogram  Exponentiation 
Francis  Kretz 

In  the  field  of  image  enhancement  simple  nonlinear  amplitude 
transformations  are  quite  useful.  For  example,  the  technique  of  histo¬ 
gram  equalization  [  1-3]  has  been  shown  to  significantly  improve  the 
detail  of  low  contrast  images  such  as  X-ray  and  Earth  Resources  Satel¬ 
lite  (ERTS)  pictures. 

Examination  of  detailed  and  well-contrasted  pictures  reveals  that 
their  histograms  are  approximately  exponential  with  a  black  level  peak. 

On  the  other  hand  poor  quality  images  usually  have  a  non-exponentially 
shaped  histogram.  For  this  reason  an  investigation  was  made  into  methods 
of  "histogram  exponentiation"  for  image  enhancement. 

Probability  Density  Function  Transformations.  Let  X  denote  a 
random  variable  whose  range  is  [  0,  l]  ;  p(x)  ,  xe[  0,  l]  be  its  probability 
density;  and  F^(x)  represent  the  distribution  function  assumed  continuous. 
Also,  let  f(  • )  be  a  continuous  monotonically  increasing  function  and 
let  Y  =  f(X)  be  an  output  random  variable  with  distribution  given  by 

F  (y)  =  Pr  {  Y  Sy}  =Pr  {  f(x)  *  y  } 

Since  f(» )  is  invertible 

F  (y)  =  P  (XS  f_1(y ) }  =  F  {f-1(y)} 
y  r  x 
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Thus,  for  a  given  input  distribution  F  (x)  and  a  desired  output  distribution 
F^(y),  it  is  possible  to  determine  the  necessary  transfer  function  f(X). 

In  the  discrete  case  the  transfer  function  can  only  be  determined  approx¬ 
imately.  If  the  transfer  function  is  chosen  such  that  F  (y)  is  uniform 
the  process  is  called  histogram  equalization. 

Histogram  Exponentiation.  In  the  histogram  exponentiation  process 
the  transfer  function  is  selected  so  that  the  output  probability  distribution 
is  of  the  form 

P  (k)  =  A  exp{-a(k-l)} 

for  k  =  1,  2 . M.  The  parameter  "a"  controls  the  shape  of  the  exponential 

function.  With  a=0,  the  result  is  histogram  equalization. 

Figure  1  shows  the  effect  of  histogram  equalization  and  exponentiation 
on  an  ERTS  picture.  The  equalized  picture  in  figure  3b  is  subjectively 
improved  compared  to  the  original.  And  the  exponentiated  pictures  exhibits 
further  improvement. 
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(b)  Eqiialization 

(M  =  31.  5  ) 


(  a)  Original 


(d)  Exponentiation 

(M  =9.4) 


(c)  Exponentiation 
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5.  Image  Data  Extraction  Projects 


Image  data  extraction  describes  the  collection  of  projects  con¬ 
cerned  with  the  detection  of  features  within  an  image  and  methods  of 
measuring  these  features. 

The  first  report  describes  an  investigation  into  a  method  for  image 
reconstruction  from  transverse-axial  density  projections  of  a  solid  object. 
The  method  utilizes  a  Fourier  transform  process  defined  on  a  polar  raster 
which  obviates  the  need  for  interpolation  in  the  transform  domain. 

A  new  project  described  in  the  following  report  is  based  upon  the 
development  of  nonlinear  optical  processing  elements.  These  elements  are 
constructed  from  halftone  transparencies  which  are  mathematically  compu¬ 
ted  and  recorded  by  a  scanning  microdensitometer. 

5.  1  Fourier- Bessel  Method  for  Transverse-Axial  Reconstruction 
Dennis  G.  McCaughey  and  Richard  P.  Kruger 

Transverse  Axial  Reconstruction  implies  the  reconstruction  of  two 
dimensional  cross  sectional  regions  of  an  object  or  signal  from  knowledge 
of  a  discrete  number  of  one  dimensional  projections.  Application  areas  in¬ 
clude  analysis  of  electron  microscope  imagery  Cl], medical  transverse  axial 
tomography  [2], and  radar  signal  analysis  [3].  Algorithms  for  this  pur¬ 
pose  may  be  divided  into  three  general  categories  :  Algebraic  [4], Con¬ 
volutional  [2], and  Fourier  Transform  domain  Cl]  methods.  The  present 
discussion  will  be  limited  to  the  latter  method. 

Transform  Processing  Techniques.  The  Fourier  method  of  reconstruction 
depends  on  the  fact  that  the  Fourier  transform  of  the  projection  is  identical 
to  the  corresponding  central  section  of  the  Fourier  transform  of  the  density 
function  C5],  When  the  Fourier  transform  operation  is  performed  by  a  digi¬ 
tal  computer,  the  transform  of  the  original  density  is  obtained  at  discrete 
points  in  frequency  space  in  polar  coordinates.  Since  few  of  these  points  will 
correspond  to  rectangular  coordinates  required  for  the  inverse  transform 
operatioi  ,  interpolation  is  required.  Also,  since  the  forward  or  inverse  dis¬ 
crete  Fourier  transform  operation  produces  the  original  function  at  this  sam¬ 
ple  point  to  within  the  truncation  errors  of  the  machine  used,*  it  would  seem 
reasonable  that  some,  but  certainly  not  all,  of  the  artifacts  present  in  the 

*  It  can  be  shown  that  the  forward-inverse  sequence  of  Fourier  transform 
operations  produce  the  original  function  exactly  at  the  sampling  points  with¬ 
out  consideration  of  the  sampling  theorem.  However,  this  is  not  to  imply 
that  the  sampling  theorem  is  without  importance. 
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reconstructed  image  are  due  to  the  particular  interpolation  process  used. 
Other  artifacts  in  the  reconstructed  cross-sectional  image  result  from 
Gibbs  phenomenon,  high  frequency  components  (due  to  edge  effects)  present 
in  the  estimated  Fourier  coefficients,  and  possible  undersampling.  While 
much  has  been  said  with  regard  to  the  sampling  theorem  concerning  image 
degradation  and  resolution  C 1 8 ,  19]  more  effort  is  necessary  to  determine 
the  relative  importance  of  the  various  interpolation  methods  in  the  area  of 
image  degradation.  Notwithstanding  the  lack  of  knowledge  concerning  the 
relative  importance  of  interpolation  as  an  important  contribution  to  image 
degradation,  it  would  seem  desirable  to  employ  an  algorithm  that  utilize 
the  Fourier  coefficients  directly  on  polar  coordinates.  Crowther  De  Rosier 
and  Klug  proposed  such  an  algorithm  using  the  Fourier- Bessel  transform 
t21  j.  Crowther  et.  al.  Cl6]  have  employed  this  algorithm  in  the  recon¬ 
struction  of  images  obtained  in  electron  microscopy  with  reasonable  re¬ 
sults.  However,  their  format  was  such  that  projections  were  not  obtained 
at  evenly  spaced  angles  which  necessitated  an  interpolation  process  to  fill 
the  polar  raster.  In  developing  an  algorithm  it  is  reasonable  to  assume  that 
projections  should  be  available  at  evenly  spaced  angles  to  avoid  this  inter¬ 
polation. 

As  a  basis  for  the  analysis  let  F(R,  9)  denote  the  Fourier  transform 
of  the  projection  at  angle  P  and  f(i',cp)  be  the  tomographic  section  in  polar 
coordinates.  The  following  then  results  [8] 

2rr  oo 

f(r,cp)=f  f  F(R,  0)exp(-2rr  j  rR  cos(P-cp)RdRd6  (1) 

0  T) 


Note  that 


cos  (9-cp)  =  sin(9-cp+rr/2). 


and 

00 

exp(ja  sin  x)  =  Jk(a)  exp 

k=  -oo 

th 

where  J^a)  is  the  k  order  bessel  function  of  the  first  kind.  By  inserting 
these  two  relationships  into  eq.  (1)  and  rearranging  the  order  of  summation 


and  integration  it  is  possible  to  obtain 

®  2tt  ® 

f(r,cp)  =  ^jkf  f  F(Rie)Jk(2nRr)exp(jkG)RdRde  exp(-jkcp)  (2) 

k=-oo  *'0  •'0 


If  f^(r)  is  defined  as 


2TT  00 


fk,r)  =  ikf  f 

Jo  Jo 


F(R,  e)Jk(2TTRr)exp(jke)RdRdP 


(3) 


it  is  possible  to  obtain  the  Fourier  series  expression 


f(r,cp)  fk(r)exp(-jl- 


cp) 


(4) 


k- 


Note  that  eq.  (4)  is  indeed  a  Fourier  series  for 


1  2tt 


fk(r)  = 


2tt 


l 


f(r,cp)exp(jkcp)dcp 


(5) 


If  constraints  to  discrete  angles  evenly  spaced  over  (0,  2tt)  are  assumed, 
eq.  (4)  becomes  a  discrete  inverse  Fourier  transform.  Knowledge  of  fk(r) 
will  then  permit  an  exact  reconstruction  of  f(r,cp)  along  radii  equally  spaced 
over  (0,  2tt).  The  factor  F(R,  8)  is  in  general  not  available  for  all  R  and  8 
since  projections  are  available  for  only  a  discrete  set  of  angles  and  because 
a  DFT  on  each  projection  produces  the  transform  only  at  discrete  fre¬ 
quencies.  Therefore  Fg(r,  8)  is  the  2-D  sampled  form  of  the  Fourier  trans¬ 
form  F(R,  8)  given  by 


Fg(R,8)  = 


2ttAR 


N 


2 ^  4^  F(iAR, 


^p)  6(R-iAR,6-  ^p) 


(6) 
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where  N  denotes  the  number  of  projections  and  I  represents  thr  number 
of  points  in  each  projection.  Inserting  eq.  (6)  into  eq.  (3)  one  obtains  an 
estimate  of  fk(r)  termed  f^r)  in  the  form  of  a  Riemann  sum 


Vr>  = 


2ttAR 


1-1  N-l 


2n  n  .  .k. 


n  £  £  F(Ri»  ~N  jk(2tt  Rir)Riexp(jk2nn/N) 
i=0  n=0 


(7) 


where  i  AR  has  been  replaced  by  R..  The  estimated  reconstruction 


is 


then 


f(r,°P)  =  ^2  f,(r)exp(-jkcf) 


(8) 


k=-°° 


Since  f(r,(#  is  a  spaced  bounded  continuous  function,  F(R,  6)  is  an  ana¬ 
lytic  function  [2]  and  thus  eq.  (7)  always  exists.  Equation  (7)  is  in 
reality  rectangular  integiation  and  if  AR  and  2tt/N  are  small,  ^(r)  is  an 
accurate  estimate  of  f^( r) .  The  ability  to  estimate  the  tomographic  section 


f(r»<#with  f(r,  cp)  is  thus  limited  only  in  the  accuracy  of  the  fk(r).  It 
be  easily  shown  that  eq.  (8)  can  be  expressed  as 


can 


00 

1-1 

N-l 

E 

J2  f(r 

* 

o 

n 

X 

i=0 

n=0 

00 

1-1 

N-l 

4-  Y*  ,k.-k  2ttAR 
^  J  e,N 

E 

EF<Ri- 

k=0  k 

i=0 

i=0 

2m 


(9) 


where 


£k  = 


2 

1 


k  =  0 
k  i  0 
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By  similar  means  it  can  be  also  shown  that 


00 


f(r,Cp)=  2Re 


ZvL  R 
ekN' 


1-  1 

£ 


N-l 

^  ^  F^Ri»  "n- )  J^(2TTR^r)R.expCjk2nn/N]exp-5k(cp  -tt/2)] 
n=0 


(10) 


where  j  is  replaced  by  exp  tjk4n/2)]4  Equation  10  can  be  computed  by  two 
fast  Fourier  transform  operations  followed  by  a  Bessel  function  weighting 
operation  and  a  final  fast  Fourier  transform  operation.  Current  results 
indicate  that  test  images  can  be  reconstructed  from  32  projections  with 
64  sample  points  in  each  projection  in  a  few  seconds  on  a  standard  comput¬ 
er.  An  example  of  a  reconstructed  disk  is  shown  in  figure  1. 

Removing  the  interpolation  process  results  in  not  only  a  significant 
reduction  in  computing  time  but  some  insight  into  a  quantitative  measure¬ 
ment  of  the  accuracy  of  the  algorithm  through  eq.  (7).  Furthermore,  re¬ 
taining  a  polar  format  produces  an  algorithm  with  the  highest  resolution 
near  the  axis  of  rotation.  Upon  first  consideration  this  may  seem  a  disad¬ 
vantage,  however,  a  higher  resolution  may  t*  obtained  in  a  region  or  interest 
simply  by  centering  the  axis  of  scan  rotation  in  that  region.  This  may  result 
in  the  need  for  fewer  projections. 

This  investigation  would  seem  to  indicate  that  the  Fourier-Bessel 
method  has  several  advantages  over  the  more  conventional  transform  method. 
It  is  still  necei  sary  to  determine  the  quality  of  the  estimate  ?k(r)  and  to  de¬ 
velop  methods  to  improve  this  estimate.  It  was  noted  earlier  that  eq.  (7)  was 
in  reality  a  rectangular  integration  formula.  A  simple  method  to  improve  the 
accuracy  of  f^( r)  would  be  to  utilize  a  trapezoidal  formula.  This  would  require 
no  increase  in  the  number  or  sample  pointf .  Also  if  the  projections  are  sam¬ 
pled  at  the  Nyquist  rate  (implicit  in  any  method),  it  can  be  shown  that  the  func¬ 
tion  can  exhibit  no  more  than  one  zero  crossing  between  the  two  adjacent  sam¬ 
ple  points.  This  would  imply  that  a  trapezoidal  formula  would  be  reasonably 
accurate.  The  algorithm  also  may  be  extended  to  images  of  much  higher  res¬ 
olution  since  fast  Fourier  transform  algorithms  are  certainly  applicable  to 
higher  dimensionality  than  currently  used.  The  main  effort  would  be  directed 
towards  calculating  the  Bessel  functions  of  higher  order.  No  particular  dif¬ 
ficulty  is  anticipated,  since  sufficiently  accurate  large  argument  approxima¬ 
tions  are  available  [9]  along  with  a  recursive  algorithm  for  small  arguments 
ClO]. 
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5.2  Nonlinear  Optical  Image  Processing  With  Halftone  Screens 
Stephen  R.  Dashiell  and  Alexander  A.Sawchuk 

Under  ordinary  circumstances  optical  image  processing  systems  • 
are  capable  of  performing  only  linear  operations  on  input  images.  There 
is  a  large  class  of  nonlinear  operations  which  could  be  of  great  value  if 
they  could  be  easily  performed  in  an  optical  system.  Among  these  are 
homomorphic  filtering  Cl],  histogram  equalization,  level  slicing,  high 
intensity  pass,  low  intensity  pass,  intensity  band  stop,  and  others  [2], 
Kato  and  Goodman  [3]  have  successfully  utilized  commercially 


available  halftone  screens  to  perform  logarithmic  and  exponential  transfer 
functions  in  an  optical  system,  thereby  allowing  them  to  filter  a  multipli¬ 
cative  noise  component  with  much  greater  success  than  if  ordinary  linear 


filtering  had  been  performed. 

If  diffraction  orders  of  the  halftoned  image  higher  than  the  on-axis 
zero-order  diffraction  component  are  used,  and  if  specially  made  halftone 
screens  are  available,  nonlinear  transfer  functions  considerably  more 
complex  than  the  monotonic  logarithm  and  exponential  can  be  obtained. 

To  understand  how  a  non-monotonic  transfer  function  might  be  ob¬ 
tainable,  consider  first  the  Fourier  transform  of  a  rectangular  array  of 
opaque  squares  of  side  b  and  center  to  center  spacing  a,  given  by 

fin  no 

2 


F(_t  (x,  y)}  =  6(fx,  fy)  - 


TT 


n=-“  m=-®  7  '  /  '  / 


(1) 


where  m  and  n  are  numbers  identifying  the  diffraction  order  impulses  in 
the  spatial  frequency  plane. 

The  transmittance  function  assumed  by  eq.  (1)  is  typical  of  what 
would  be  obtained  if  a  constant  transmittance  were  halftoned  and  photo¬ 
copied  on  high  contrast  copy  film.  The  presence  of  the  sine  terms  indicates 
that  non-monotonic  behavior  could  be  expected.  This  becomes  more  evident 
in  the  special  case  where  a  (0,  1)  order  specified  by  m  =  1,  n  =  0  is  selected. 
In  this  case  eq.  (1)  reduces  to 


F{t(x,y)](o(i)  = 


(2) 


Inverse  Fourier  transforming  this  expression  and  squaring  to  get  the  out¬ 
put  intensity  yields 


(3) 


as  the  final  intensity  output. 
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This  result  is  useful  because  the  value  of  b  is  a  function  of 
the  input  intensity  to  the  system,  or  equivalently,  is  a  function  of  the 
density  of  an  input  transparency.  In  operation,  the  Fourier  transforming 
and  inverse  transforming  are  instantaneously  performed  in  a  coherent  op¬ 
tical  parallel  processor  [5],  and  selection  of  the  diffraction  orders  is  ac¬ 
complished  by  simple  spatial  filters.  These  filters  are  chosen  to  pass  the 
low  spatial  frequency  information  in  the  original,  thus  desampling  the  half- 
toned  picture  while  the  nonlinear  operation  is  performed. 

To  illustrate  how  an  operation  such  as  level  slicing  could  be  per¬ 
formed  using  this  first  diffraction  order,  a  halftone  screen  wa  3  made  by 
photographing  crossed  Ronchi  rulings.  The  resulting  screen  had  trans¬ 
parent  squares  on  a  partially  transmitting  background.  The  sides  of  the 
squares  were  one-half  the  distance  between  squares.  The  subject  trans¬ 
parency  was  then  photographed  through  this  halftone  screen,  and  three  dif¬ 
ferent  regions  were  then  present  on  the  copy  film  after  development.  In 
the  first  region,  where  the  subject  was  sufficiently  dense  that  the  copy  film 
did  not  expose  even  through  the  clear  squares,  there  were  no  dots  or  b  =  0. 
In  the  second,  where  the  subject  was  less  dense  such  that  the  copy  film  did 
expose  through  the  clear  squares,  but  not  through  the  darker  background, 
the  dot  parameters  were  b  =  a/2.  Third,  where  the  subject  was  still  less 
dense  such  that  the  copy  film  exposed  through  the  clear  squares  and  through 
the  background,  the  parameters  were  b  =  a. 

The  cases  b  =  0  and  b  =  a  give  IQUT  =  0,  however  b  =  a/2  yields 

*OUT  a  /n  ^0.  Thus,  only  the  range  of  densities  which  gave  dots  appear 
in  the  output  and  densities  above  or  below  this  range  do  not  appear.  This  has 
beer  experimentally  verified  by  slicing  a  picture  into  6  different  levels,  us¬ 
ing  a  single  halftone  screen  but  varying  the  exposure  time  on  the  copy  film. 
The  position  of  the  sliced  level  is  controllable  by  varying  exposure  while 
the  width  of  the  level  is  fixed  by  the  density  difference  on  the  screen  be  - 
tween  the  clear  squares  and  the  darker  background.  This  technique  for  ob¬ 
taining  nonlinear  transfer  functions  is  not  limited  to  level-slicing  and  other 
simple  functions.  Using  digitally  produced  halftone  screens  with  appropri¬ 
ately  selected  density  profiles,  more  complex  operations  are  possible. 

Although  these  preliminary  experiments  in  nonlinear  processing  in¬ 
volve  photographic  processing  and  coherent  optical  techniques, they  are  es¬ 
pecially  well  suited  to  applications  in  hybrid  digital/optical  syst2ms  [4]. 
Several  recl-time  optical  input  devices  with  adjustable  parameters  and 
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thresholding  are  now  under  development,  and  combing  these  with  the  half¬ 
toning  operation  would  make  possible  a  real-time  nonlinear  optical  parallel 
processor.  Such  a  system  would  be  relatively  inexpensive  and  simple,  and 
would  avoid  problems  of  scanning  and  display.  Electrically  controllable 
selection  of  diffraction  orders  in  the  Fourier  plane  of  the  system  would  per¬ 
mit  fast  modification  of  the  transfer  functions.  Immediate  future  work  will 
involve  the  production  of  several  types  of  special  halftone  screens  in  the 
IPI  plotting  microdensitometer. 


References 


Oppenheim,  R.  W.  Schaefer  and  T.  G.  Stockham,  Jr.  ,  "Non¬ 
linear  Filtering  of  Multiplied  and  Convolved  Signals,"  Proc.  IEEE 
Vol.  56,  No.  8,  pp.  1264-1291,  August,  1968.  - 1 

H.  C.  Andrews,  A.  G.  Tescher  and  R.  P.  Kruger,  "Image  Proc¬ 
essing  by  Digital  Computer,"  IEEE  Spectrum.  Vol.  9.  No.  7 
pp.  20-32,  July,  1972.  - C - 1 

H.  Kato  and  J.  W.  Goodman,  "Nonlinear  T ran? formations  and 
Logarithmic  Filtering  in  Coherent  Optical  Systems,"  Opt.  Comm. 
Vol.  8,  No.  4,  pp.  378-381,  August,  1973.  - 


4. 


5. 


GUeS‘  Edit°r'  special  is8ue  on  Optical  and  Digital 
KUy™<S!e° Dm"8’"’8’  gg.ical  En8ineerinK,  Vol.  13,  No.  3, 


J.  W.  Goodman, 
1968. 


Intro,  to  Fourier  Optics. 


McGraw-Hill,  New  York, 


-68- 


The  image  analysis  projects  are  concerned  with  the  back¬ 
ground  technology  necessary  to  effectively  design  image  coding,  res¬ 
toration,  enhancement,  and  data  extraction  systems.  Of  particular 
interest  are  models  of  the  human  visual  system  for  monochrome  and 
color  images,  and  the  development  of  quantitative  measures  of  image 
fidelity  and  intelligibility. 

In  the  first  project  the  continued  research  on  a  model  of  human 
color  vision  is  described.  This  model  is  nearly  completely  developed 
and  tested,  with  success.  The  first  applications  of  the  model  are  now 
int  roduced. 

The  next  report  considers  the  development  of  a  more  accurate 
model  for  human  monochrome  vision.  This  model  incorporates  a  linear 
filtering  element  before  the  photoreceptors  to  account  for  the  optical 
degradations  of  the  eye. 

6.  1  A  Quantitative  Model  of  Color  Vision 
Werner  F rei 

The  first  part  of  this  report  summarizes  recent  results  of  an  effort 
to  describe  some  major  properties  of  human  color  vision  with  a  simple 
neuro-physiological  model.  The  second  part  investigates  an  application 
of  the  model  to  the  optimal  quantization  of  color  image  signals. 

Neuro-physiological  model.  The  visual  model  in  question  is  the 
result  of  a  broad  study  of  the  neuro-anatomy  of  the  human  eye  and 
psychophysics  of  vision  C  1,  Z]  .  Figure  1  shows  a  block  diagram  of  the 
model  which  consists  of  a  layer  of  photoreceptors  connected  to  a  layer 
of  summation  cells.  The  photoreceptors  comprise  three  distinct  classes 
corresponding  to  Judd's  fundamental  sensations,  with  a  quasi-logarithmic 
response  to  light  stimulation.  Three  distinct  patterns  of  connection  are 
assumed  between  the  two  layers,  such  that  the  output  of  each  summation 
cell  corresponds  to  one  of  three  opponent  sensations  "dark-light"  "red- 
green"  or  "  blue-yellow.  " 


-69- 


CL 

<D  CD 
U  ^  U 
O  rj 
<D  — 
o  O  Cl 
ru  w 

II 


<L>  Q.  ZJ 

O  a,  Q. 
O  Q  *•“ 
CLOU  3 
CO  i-  O 

II 


Is 

(/)  o 

C  CL 

+-  CO 
II 

Vx 


O  CD 

cO'u 
-  o  o 

o>  E  Q- 

=  co 


2^o 

o  4-  o 

JC 

Z  O  CL 


I 

o 

C  O 

P  o  JO 
C  •—  r> 


§J£s 

O  o  .£  2 


CL  (S> 
CD 


-:o- 


At  high  intensities  and  moderate  saturations,  the  receptor 
response  functions  may  be  simplified  to 

t.  =c.  In  (t. )  t  >  0  1 5) 

1111 

for  t.  >0.  The  spaces  of  color  tri  stimulus  7  and  perceptual  quantities 
^  >  respectively,  are  defined  to  be  vector  spaces  with  algebraic  laws 
of  composition  defined  as  follows 


Ta-LTB  (tlAtlB,t2AtlB,t3At3B) 


T  =  <vy;>T 


GA—  GB  <glA+glB’ g2A+g2B’ g3A+g3B) 


r  Jl_  G  =(rg^,  rgz>  rg^ 


( 6a) 

'6b) 

(6c) 

(6d) 


where  T^.T^,  T  are  *he  tristimulus  vectors  of  arbitrary  colors  and  G  , 

A 

GB’G  are  the  corresPonding  neural  quantities.  The  symbol  _]_  denotes 
vector  addition  and  _IL  multiplication  by  scalars  r  (real  numbers).  The 
nonlinear  transformation  defined  by  the  model  can  now  be  easily 

verified  to  be  a  linear  vector  mapping  of  the  above  vector  space  T  into 
by 


H©  N(T  !  T  )  =  H©NT  +  H®  NT 

A  —  B  A  B 

(7a) 

H©  N(r  JJ_  T)  =  rH©  NT. 

(7b) 

This  definesa  generalized  law  of  superposition  of  colored  sensations. 
For  example,  multiplying  (  II  )  the  tristimulus  values  of  a  color  by 
a  scalar  r  aj  defined  by  eq.  (6b)  yields  a  color  appearing  approximately 
r-times  as  bright,  sa  urated,  etc. 

The  analysis  of  the  visual  response  to  an  arbitrary  two-dimen  - 
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(2c) 


t  =  JU  S(X)  t..  dX 

i  ,  lX 

\ 

where  forj=  1,  2,  3  are  independent  linear  transformations  of  the 

CIE  Xy  '  yX  '  7X  C°loT  matchin8  functions;  X,  and  X  are  the  limits  of 
the  visible  spectrum. 

The  mapping  N  describes  the  neural  response  of  the  receptors, 
assumed  to  be  a  quasi- logarithmic  function  of  the  respective  light 
energies  absorbed 


t.  =  c.  In  fk.t.  +t  ) 

ii  ii  iO  <3) 

where  c.,  k.,  t.Q  are  constants.  Finally  the  transformation  H(x,  y) 
represents  linear  weighted  summations  of  the  receptor  outputs  in  the 
lateral  inhibition  process. 

For  simplicity  three  distinct  shift  invariant  patterns  of  summations 
are  assumed,  which  produce  opponent  neural  signals  g,<x.  y),  g^x.  y)  g  (x.  y) 

corresponding  to  "dark-light. "  "red-green"  and  "yellow-blue"  perceptual 
quantities 


G(x,y)  _  Cg1(x.y),g2(x,y)g3(x,y)] 


+  00 


gl(x,y)  f  h1(x-?,y-ri)c1ln(t1(5.y)+  t^Jd^q 


-.00 
+  00 


82<x.y) 


d  rdr| 


-.00 

+00 


83.x.y,  .//h 3(,-e  .MIVn(ffj 


d§dr) 


(4a) 


(4b) 


(4c) 


(4d) 
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The  three  opponent  sensations  form  the  basis  of  a  "perceptual" 
euclidian  color  space.  In  that  space  "brightness"  is  equal  to  the  norm 
of  the  vector  sum  of  the  three  opponent  signals. 

The  pattern  of  connections  between  the  two  layers  also  predicts 
spatial  contrast  phenomena.  For  intensities  and  color  saturations  of 
practical  interest,  tne  spaces  of  color  tristimulus  values  and  neural  signals, 
respectively,  are  defined  as  algebraic  vector  spaces  with  distinct  laws  of 
composition  and  a  linear  vector  mapping  from  the  first  into  the  second. 

From  this,  a  generalized  law  of  superposition  of  colored  sensations 
is  defined  which  enables  fast  computational  procedures  to  predict  the 
perception  of  color  in  complex  visual  fields.  The  mathematical  structure 
of  the  model  is  now  briefly  outlined. 

The  model  defines  a  mapping  of  the  retinal  image  space  .'into  a 
space  of  opponent  neural  signals  If,  with  SG  ,  x,  y)  e  J>  de  signaling  the 
spectral  energy  distribution  of  light  striking  the  retina  and  Gfx,  v)  e 
a  vector  valued  "perceptual"  quantity;  x.yare  geometric  coordinates 
in  the  retinal  image  plane.  The  mapping  consists  of  a  cascade  of  a 
linear  (T^),  a  nonlinear  (N)  and  a  second  linear  <H(x,y))  transformation 
as  given  by 


G(x.y)  =  H(x,y)©N  T^SfX.x.y) 

where©  denotes  the  convolution  integral.  Whenever  possible,  the 

geometric  coordinates  (x,y)  will  be  omitted  to  simplify  the  notation. 

The  linear  transformation  T^  maps  the  spectral  energy  distribution  of 

light  S(X)  into  a  colorimetric  tristimulus  space  T .  The  factor  T  des- 

\ 

cribes  the  quantum  absorption  mechanism  in  three  types  of  photopic 
receptors  and  is  defined  for  the  standard  observer  as 


T^:  S(X,  x,  y)-*T(x,  y) 
T  =(‘1*‘z-‘3>T 


(2a) 

(2b) 
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Sion**  field  can  now  be  cast,  in  a  first  approximation,  into  the  frame¬ 
work  of  linear  systems  analysis.  In  particular,  spatial  contrast  ef¬ 
fects  and  resolution  are  predictable  by  a  linear  two-dimensional  vector 
filtering  V(ju,  jv)  ini'-  space 


G(x,y)  =  3'  {V(ju,  jv)  <7[  CNT  ^S(X,  x,  y) ]  } 


(8) 


where  .7 and  3  ‘denote  the  two-dimensional  Fourier  transformation 
and  its  inverse  respectively,  (u,  v)  are  spatial  frequencies  in  the  (x,y  ) 

direc  ion  and  V  (ju,  jv)  .s  the  Fourier  transformation  of  the  impulse  response 
matrix  H(x,  y).  The  model  has  been  previously  shown  [2]  to  be  in  agreement 
with  the  following  psychovisual  factors:  a)  metameric  color  matching,  b)  ad- 
ditivity  of  luminances,  c)  adaptation  and  simultaneous  contrast.  An  additional 

hypothesis  allows  to  predict  the  color  discrimination  ability  of  the  visual 
system. 


Quantitative  Color  Discrimination.  The  statistical  errors  of 
color  matching  experiments  can  be  approximately  predicted  by  assuming 
that  the  neural  signals  gj.g2.g3  are  corrupted  by  independent.additi ve 
gaussian  random  fluctuations. 

Extensive  studies  of  color  matching  errors  [  3]  have  shown  that 

the  loci  of  the  experimental  standard  deviations  around  fixed  reference 

colors  are  ellipsoids  in  the  CIE  x,  y,  l-  space,  with  ta  logarithmic 

function  of  the  luminance  Y  ;  {,  =  0.  2  Ioe  (Y) 

10  „ 

Let  the  probability  of  confusing  a  color  T  =rT£{\)  with  some 
reference  color  T  =  T^S(X)  be  equal  to 


p(G  I  T)  =  ~  exp  {  -y  [  G-G)T(G-G)]  1 
with  G  =  CNT,  G  =  CNT. 

Under  this  hypothesis,  the  experimental  ellipsoids  should  cor¬ 
respond  to  unit  diameter  spheres  in  -  space.  Since  the  model  has  a 
separate  luminance  channel,  luminance  and  chromaticity  errors  can 
be  considered  separately.  (Brown  and  MacAdam's  data  [  3]  suggests 
generally  small  correlations  of  luminance  and  chromaticit/  errors.) 


(9) 
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Using  the  simplified  expressions  for  and  g^  at  high  radiances, 
and  letting  =  k^  =k^  =1  (see  discussion  on  adaption  below),  the  semi¬ 
axes  of  MacAdam's  25  ellipses  were  mapped  into  J'-  space,  adjusting 
the  constants  c^  and  c^  to  obtain  a  best  fit  with  unit  diameter  circles 
in  a  plane  g^  =constant.  Table  1  shows  the  lengths  of  the  mapped 
semi-axes  ini'-  space  for  c^  =123.2  and  c3  =18.8.  Figure  2b  shows  the 
result  of  the  inverse  operation:  unit  diameter  circles  ini'-  space, 
gj  =  constant  ,  mapped  into  the  CIE  x-y  chromaticity  diagram.  The 
constant  c^  =67.4  is  determined  from  the  average  standard  deviations 
of  luminances  reported  in  the  same  study  [  3] 

cf  A  t{  j  2  =  60.  104  A  yj  2  «  60. 10  W  (10) 

Color  Distance.  The  statistical  deviations  of  color  matches 
provides  a  measure  of  the  visual  threshold  color  discrimination.  Practice 
shows  that  Just  Noticeable  Differences  (JND's)  are  approximately  equal 
to  three  times  the  standard  deviations  [  4*1  . 

According  to  Schrodinger' s  hypothesis  [  4]  the  perceptual  differ¬ 
ence  |  A  B  |  between  two  arbitrary  colors  A  and  C  can  b<  determined 
as  the  smallest  number  of  JND's  separating  the  two  colors 

C 

|4B|=(t/d) 

where  dS  is  a  line  element  and  the  path  cf  integration  is  a  geodesic  line 
[4]. 

Under  the  assumption  of  neural  channels  perturbed  by  independent 
random  fluctuations,  a  line  element  can  be  derived  from  the  present 
model 
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where  one  tristim- 


«I=[E  (E-^d,.)*]* 


(\2a) 
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and,  at  high  radia  nee s,  approximately 


111 
-c.  +  c  4c  _ 


'  “ 1  '^T  ‘(vT  ‘ 


klk2tit2  dt*dt2 


, .  .  dt  dt 

WlS  1  3- 


Two  particular  cases  are  of  interest,  namely  the  line  elements  for 
constant  luminance  and  achromatic  differences  respectively 


Y--  Cst 


c  dt 

I  =-! — L 

achromatic 


(I2e) 


1  quation  (12e)  shows  that  the  model  is  consistent  with  Weber’s  law. 

Examination  of  eq.  (12a)  reveals  that  the  geodesics  are  straight 
lines  in>-  space.  The  measure  of  perceptual  color  d  fference  thus 
reduces  to  a  simple  euclidian  metric 


-78- 


(13) 


C 

I4BI  .-L/dS.-L  [<cA-  cc,T(GA-Gc)]^ 

A 

The  >-  space  can  therefore  be  called  an  approximately  euclidean  color 
representation  space. 

Brightness,  Lines  of  Constant  Hue  and  Saturation.  Assume  that 
the  rensation  of  brightness  produced  by  a  color  A  is  proportional  to 
the  perceptual  difference  between  that  color  and  blacK  Brightness  is 
then  proportional  to  the  norm  of  G 


B 


M4) 


The  failure  of  brightness  additivity  in  color  mixtures  then  derives 
immediately  from  Schwartz's  inequality 


B 


A+C 


=  'G.*V  1*8  *Br  -IO.I+  |G  I 


(15) 


The  model  therefore  predicts  that  the  brightness  of  color  mixtures  Is 
always  smaller  or  equal  to  the  sum  of  the  component's  brightnesses. 

Schrodinger' s  hypothesis  [  4]  that  colors  of  constant  saturation 
lie  on  geodesic  circles  with  constant  luminance  will  now  be  explored. 

In  A  -  space,  such  geodesics  are  ordinary  circles  in  a  plane  g^  =  constant 
and  centered  at  g^  =g^  =  0.  The  shortest  distanc?  from  any  point  on 
one  circle  to  a  point  on  another  circle  is,  of  course,  along  a  straight 
line  pas  sing  through  g^  =%2  ~  Radial  lines  in  -  space  would  then, 
according  to  Schrodinger,  be  projections  of  lines  of  constant  hue  on  the 
plane  g^  =  constant.  Figure  3b  shows  a  set  of  equidistant  concentric 
(irclesand  radial  lines  in  A-  space  mapped  into  the  CIE  x-y  diagram. 
These  lines  resemble  very  much  the  lines  of  constant  saturation,  and 
especially  constant  hue  found  in  empirical  color  classifications  (compare, 
for  example,  with  the  Munsell  system  [  l]  ).  It  seems  therefore  that 
the  perception  if  hue  is  governed  by  the  ratio  of  the  chromatic  signals 
^2  ar>d  g-^  and  the  sensation  of  color  saturation  would  derive  from  the 
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figure  6.1-3.  Lines  of  constant  hue  and  constant  saturation 
(Schrodinger's  hypothesis) 

80 


J 


T  - 

norm  of  a  vector  sum  ^0,  g^f  g^)  (0,  g^,  g^)]  2  of  the  same  signals. 

Applications.  It  has  been  shown  r  5]  that  the  use  of  a  simple  model 
of  achromatic  vision  can  improve  considerably  the  efficiency  of  coding,  en¬ 
hancement  and  restoration  techniques  of  black  and  white  images.  This  is 
essentially  done  as  follows:  let  I(x,  y)  represent  the  intensity  of  an  image 
and  M  a  mapping  into  a  perceptual  space  (here  G  =  H  *  NI(x,  y)  is  a  scalar 
function  of  I).  Image  processing  takes  place  in  1  -  space  and  an  inverse 
mapping  M  is  performed  before  displaying  the  processed  imap;. 


This  idea  has  been  extended  to  color  image  processing  for  a 
problem  only  apparently  trivial,  namely  the  efficient  quantization  of 
color  image  signals  [  6,  7]  .  A  color  image  has  been  mapped  into 
■£>  -  space  according  to  eqs.  (  3  )and(4)with  k.  =1,  t  =0,  h.(x, y)  =1, 
i=  1,  2,  3.  The  signals  g  ,  g  ,  g  were  then  quantized  linearly  such  that  each 
quantization  step  corresponds  to  approximately  three  times  just  notice¬ 
able  differences  iJND).  Accordingly,  g^  was  quantized  with  32  levels, 
g  fourteen  levels  and  g  twelve  levels.  The  quantized  signals  were 
then  transformed  back  into  red,  green  and  blue  signals  and  displayed  on 
a  color  TV  monitor.  Figure  4a  and  b  show  the  original  and  quantized 
picture  respectively.  Note  that  artificial  contouring  is  almost  inexistant 
despite  the  coarse  quantization  used.  Figure  4c  and  4d  show  the 
subjective  effect  of  linear  quantization  on  the  signals  g^  and  g^ 

'gj  ^constant).  The  boundaries  of  the  color  diagram  is  determined  by 
the  primaries  of  the  display  device,  and  subjective  color  variations 
across  the  diagram  are  almost  uniform  (except  in  the  neighborhood  of 
the  red  primary  which  lies  very  close  to  the  spectrum  locus).  The 
quantization  of  the  chromatic  component,  in  figure  4b  corresponds  to  that 
of  figure  4d. 
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b.2  A  New  Look  at  the  V;sual  System  Model 
Hideo  Murakami  and  Ernest  L.  Hall 

Study  of  the  human  visual  system  and  visual  perception  is 
important  in  many  fields  especially  psychology,  anatomy,  and  ph/siology; 
however,  a  knowledge  of  visual  perception  is  a  fundamental  requirement 
for  the  designer  of  picture  processing  techniques.  In  fact,  the  basic 
difference  between  the  numerical  analysis  of  matrices  and  computer 
picture  processing  is  based  on  the  psychovisual  concept  of  a  visual  image. 

A  new  visual  system  model  i?  presented  here.  The  structure 
of  the  model  is  based  on  the  histological  and  anatomical  structure  of  Lie 
human  visual  system.  A  detailed  description  of  the  human  visual  system 
may  be  found  in  [  1-3]  .  Furthermore,  the  analogous  optical  system  for 
the  occular  media  is  well  documented  ("  4-c]  .  The  currently  accepted 
spatio-temporal  model  is  briefly  presented  and  was  constructed  from  a 
review  of  current  literature.  Then  the  proposed  visual  system  model 
is  presented  and  the  spatial  and  temporal  responses  are  explored  in 
detail.  An  analysis  of  the  combined  spatio-temporal  response  and  a 
conjecture  about  color  response  are  also  described. 

Present  Spatio-  Temporal  Model.  An  overall  model  of  the  visual 
system  was  not  found  in  the  literature  although  several  references 
emphasizing  the  spatial  and  temporal  frequency  characteristics  of  the 
human  visual  system  are  available.  From  a  careful  survey  of  these,  one 
may  arr  ve  at  what  may  be  called  the  "presently  accepted  visual  system 
model.  The  spatial  elements  of  this  model  consist  of  a  nonlinearity 
followed  by  a  bandpass  spatial  filter.  It  should  also  be  noted  that  if 
one  assumes  a  retinal  stimulus  image  as  does  Stockham  [  7]  ,  the  above 
elements  do  not  conflict  severely  with  the  model  proposed  in  this  paper. 
However,  in  experimental  work  the  stimulus  image  is  presented  to  the 
occular  media  not  to  the  retina.  The  shape  of  the  nonlinearity  is 
commonly  taken  as  logarithmic  although  other  forms  are  also  found. 

The  shape  of  the  bandpass  characteristics  varies  from  strictly  low  pass 
to  bandpass  but  generally  a  bandpass  is  used.  The  temporal  model  as 
well  as  the  combined  spatio-temporal  model  is  usually  considered  as  a 
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nonlinearity  followed  by  bandpass  temporal  or  spatio-temporal  filters. 
This  presently  accepted  model  or  subsets  of  it  are  described  by  Corn- 
sweet  [  l]  ,  Davidson  f  6]  ,  Robson  [  8]  ,  Yasuda  and  Hiwatasi  [  10]  , 
etc.  This  basic  model  has  been  found  useful  in  several  experiments; 
however,  certain  experimental  difficulties  especially  at  high  spatial 
frequencies  have  arisen,  e.  g.  ,  Davidson,  Cornsweet. 

The  authors  propose  that  ignoring  the  occular  media  of  the  eye, 
as  is  done  in  this  representation,  produces  a  fundamental  flaw  in  this 
basic  model  which  has  l^d  to  the  experimental  difficulties.  Further¬ 
more,  it  will  be  shown  that  a  new  model  can  be  developed  which  is  no 
more  complex  than  the  present  model,  is  based  on  the  histological 
and  anatomical  structure  of  the  visual  system,  predicts  the  physiological 
function  of  the  visual  system,  and  is  useful  in  predicting  psychophysical 
phenomena  as  shown  by  calculations  and  comparisons  with  well-known 
psychophysical  measurements. 

Proposed  Visual  g,,stem  Model.  Perception  takes  place  in  the 
brain.  This  non-obvious  fact  has  been  proven  by  observing  persons  who 
completely  lost  their  sight  from  damage  to  the  visual  cortex  without 
damage  to  any  other  part  of  the  visual  system.  Therefore,  a  complete 
model  of  the  visual  system  must  include  models  of  the  eye,  the  optic 
tract,  and  the  visual  cortex.  Although  the  state  of  development  of  a 
complete  visual  system  model  is  still  in  its  embryonic  stages,  the 
authors  believe,  as  Cornsweet  states,  that  "it  is  pleasing  to  consider  the 
possiblity  that  the  human  nervous  system,  despite  its  obvious  overall 
complexity,  is  really  composed  of  a  very  large  number  of  repetitions 
and  slight  variations  of  a  few  simple  mechanisms"  [  l"|. 

A  simplified  functional  description  of  the  perceptual  process  may 
be  given  as  follows.  The  stimulus  imige  is  generated  from  a  scene  and 
transmitted  to  the  image  forming  elements  of  the  eyes.  These  elements 
focus  the  image  on  the  retinal  receptors.  The  retinal  neural  connections 
may  then  perform  certain  processing  operations  on  the  neural  signals. 
The  next  element  is  the  neural  processing  performed  in  the  optic  neural 
connections  from  the  eyes.  Cortical  processing  tasks  include  memory 
reference  images,  various  types  of  information  processing,  features 
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extraction,  decisions,  or  descriptions  of  various  elements  of  the  scene. 
The  results  of  this  process  may  be  indicated  in  several  ways  by  the 
response.  The  complexity  of  the  perception  process,  coupled  with  the 
desire  to  predict  even  limited  responses,  motivate  the  development  of 
simplified  models  of  the  visual  system. 

The  proposed  model  of  the  visual  system  is  shown  in  block  diagram 
form  in  figure  1.  The  stimulus  image  is  represented  by  an  image 
function.  The  image  forming  elements  of  the  eye  are  represented  by  ar 
optical  filter.  The  retinal  receptor  response  is  divided  into  two  pathways, 
one  for  rods  and  one  for  cones.  Each  pathway  contains  a  selective 
spectral  filter,  a  low  pass  temporal  filter,  and  a  low  pass  spatial  filter. 
These  are  followed  by  a  threshold  function,  a  log  conversion,  and  a  high 
pass  spatio-temporal  filter.  The  neural  connections  in  each  channel  are 
represented  by  a  spatial  optical  filter  and  summations.  The  neural 
processing  in  the  optic  nerve  is  again  represented  by  a  summation.  The 
transmission  link  to  the  cortex  may  be  represented  by  a  set  of  transmission 
lines.  Finally,  the  cortical  computations  may  be  modeled  as  summations. 
The  stage  of  development  of  the  model  of  the  visual  system  is  far  from 
complete,  although  a  large  amount  of  work  has  been  accomplished  ri-12]  . 

At  present,  certain  elements  such  as  the  im-ge  forming  mechanisms 
can  be  modeled  very  accurately;  other  elements,  for  example  the  cortical 
computations,  are  onl>  modeled  in  a  conjectural  manner.  Therefore, 
one  may  expect  to  find  insight,  but  not  a  complete  solution  from  the 
following  description  of  the  visual  model. 

Spatial  Frequency  Response.  A  fundamental  difference  between 
the  presently  accepted  model'  and  the  proposed  model  is  the  location  of 
the  elements  which  affect  the  spatial  frequency  response.  The  currently  ac¬ 
cepted  model  is  usually  represented  by  a  nonlinear  response  such  as  a 
logarithmetic  response  followed  by  a  spatial  bandpass  filter.  In  the  new 
model  the  first  element  encountered  is  a  spatial  low  pass  filter,  then  the 
nonlinearity  and  finally  a  spatial  high  pass  filter.  The  consequences  of 
these  different  configurations  will  now  be  explored  in  detail. 

The  spatial  low  pass  effect  before  the  logarithmic  response  in  the 
new  model  comes  from  the  optical  aberration  caused  by  the  lens  in  the  eye 
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Figure  6.2-1.  Backward  inhibition  model  for  two  receptors. 


and  the  light  scatter  within  the  retinal  tissue.  The  measurement  of 
this  optical  performance  of  the  human  eye  has  been  investigated  by  several 
researchers  [  4-12]  .  Despite  the  lack  of  direct  access  to  its  image 
plane  at  the  retina,  Gerald  Westheimer  and  Fergus  W.  Campbell  [  4] 
measured  the  light  distribution  formed  by  the  living  human  eye  and 
found  that  the  line  spread  function  cf  the  human  eye  is  approximately  of  the 
form  exp  1  -a  |x!  }  where  a  depends  mainly  upon  the  pupil  diameter.  For 
a  white  light  and  pupil  diameter  ol  3mm, a  =  0.  7.  From  this  line  spread 
function  one  may  easily  compute  the  modulation  transfer  function.  For 
simplicity  one  may  assume  the  one  dimensional  case  and  calculate  the 
Fourier  transform  H^(w),  to  determine  the  frequency  response  of  the 
low  pass  filter 


Hj(  'll ) 


f 


f(x)exp{  -j  tix]  dx 


which  results  in 


H^ui) 


2a 


2  ,  ,2 
a  +(ui) 


One  may  see  that  the  3dB  down  cutoff  frequency  of  this  low  pass  filter 
is  approximately  6.6  cycles/degree.  Thus,  it  can  hardly  be  disputed 
that  anatomically,  the  low  pass  filter  modifies  all  images  information 
before  the  receptors  are  encountered.  However,  other  questions  naturally 
arise.  Is  the  log-bandpass  filter  equivalent  to  the  low  pass-log-high  pass 
combination9  Is  the  low  pass  effect  of  the  occular  media  of  secondary 
importance  to  the  low  pass  effect  of  the  discrete  receptor  array  or  the 
low  pass  effect  of  mutual  receptor  excitation9  To  answer  these  questions 
the  model  of  the  retinal  neural  networks  must  be  considered. 

Retinal  Neural  Networks  .  G.  G.  Fur-nan  [  11  ]  proposed  four  types 
of  neural  networks:  a  forward  inhibition  model,  a  backward  inhibition 
model,  a  forward  shunting  model,  and  a  backward  shunting  model.  Four 
of  these  essentially  have  the  same  effect,  i.  e.  a  spatial  hi  ;h  pass  filter 
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effect.  The  backward  inhibition  model  is  employed  in  our  model  and 
is  explained  in  detail  in  [l]  . 

The  backward  inhibition  model  for  two  receptors  is  shown  in 
figure  2.  In  the  figure,  the  receptors  have  a  logarithmic  response 
to  the  incoming  light  intensity;  thus  they  correspond  to  the  logarithmic 
operation  in  the  proposed  model.  In  figure  3,  e.  is  the  frequency 
at  which  receptor  i  would  produce  pulses  if  the  receptor  i  alone  were 
illuminated  (the  level  of  excitation  of  the  receptor  i);  g.  is  the  frequency 


of  the  pulse  of  the  receptor  i  after  the  network  (pulses/sec.);  b  is 

ij 

the  inhibitory  coefficient  representing  the  strength  of  the  inhibition 


that  g,  exerts  on  g..  If  one  neglects  the  inhibitory  threshold  and  the 


case  when  frequencies  of  pulses  are  below  inhibitory  threshold,  the 
following  equations  hold 


E  -  e  -  b  g 
B1  1  12e2 


82  =e2'b12gl- 


The  general  equation  considering  other  receptors  is 

n 

g,  -  e,  “S  ' h  g  ,  i=l,2,...,n. 


g.  =  e.  -y  b.  .g. 
i  i  Z— i  ij  J 

j  =1 


The  equations  may  also  be  written  in  concise  matrix  form  with  the 
definitions 
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rzzi  inhibition  model  for  two  receptors. 
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part  I  part  2 


(a)  Proposed  model  of  spatial  characteristics 


part  I  part  2 

r - 1  ' - 1 


(h)  Modified  model  of  spatial  characteristics. 


Figure  6.2-3.  Proposed  model  of  spatial  characteristics. 
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With  this  matrix  notation  the  above  equation  will  be: 


g  =  e  -  Bg_ 


Then 


{ I  +  B)  g  =  e 


or 

g  =  (I  +  B)"1  e 


Now  assume  that  there  is  no  self-inhibitory  interaction  and  that  the  in¬ 
hibitory  interaction  is  an  exponentially  decreasing  function  of  the  distance 
of  the  receptors 


b.. 
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aQexp{  -a  [  i - j  I  } , 


=  j 


Under  this  assumption 
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This  matrix  is  Toplitz  since  all  diagonal  terms  are  equal,  and  fur¬ 
thermore  may  be  interpreted  as  corresponding  to  a  linear,  shift  invari¬ 
ant  system.  Thus,  this  convolutional  matrix  has  impulse  response 
representation 


g(x)  =  (1  - aQ )  M*)  +  aQexp(-a|  x|  ) 
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The  Fourier  transform  of  the  impulse  response  is 


G(id)  =  2r  (l-aQ)  + 


2a 

2  2 

a  +  tin 


Thus  (I  +  B) 


has  the  frequency  characteristic 


H2(UJ) 
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2  2 

2a^a+2n  (l-a^)(a  +W  ) 


This  function  represents  a  high  pass  filter  with  H  (0)  = - - - 

2  2a  +2na(l-a  ) 


1 


and  lim  H  (w)  =  - 

.  _  2  2n(l-a„) 

vr* °°  0 


From  comparing  this  function  to  the  experimental  results  of 
Davidson  [  6l  !t  appears  that  suitable  parameter  values  to  use  for  this 
comparison  are  :  a  =0.01  and  a^  =  0.6. 

Davidson  [  6]  measured  approximate  modulation  transfer  functions 
of  the  human  vision  at  several  different  levels  of  contrast.  He  called 
these  measured  functions  "describing  functions"  to  emphasize  the  nonlinear 
aspects  of  the  visual  system.  His  experiments  showed  that  the  describing 
functions  vary  with  contrast  at  high  frequencies.  This  variation  has 
been  attributed  to  the  experimental  apparatus.  It  has  been  found 
possible  to  explain  this  variation  by  a  simple  analysis  of  the  proposed 
model. 


Based  on  the  hypothesis  that  nonlinear  characteristics  of  the  human 
visual  system  could  be  traced  to  a  simple  logarithmic  transformation  at 
the  earlier  stage,  Davidson  used  the  stimuli  of  the  form, 

f(x,y)  =  f  exp{  A  cos(Uty)]  to  present  to  the  eye  in  his  experiments.  Where  f° 

is  a  constant;  A  is  a  constant  to  measure  the  contrast  of  the  signal; 

^  /2tt  is  the  spatial  frequency  of  the  signal.  The  logarithm  of  the  stimuli 
o 

will  be  log  f  +  A  cos(<i*y),  which  has  a  desirable  form  to  find  the  fre- 
quencv  characteristics  of  the  linear  filter  after  the  logarithmic  trans¬ 
formation. 
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The  experiment  was  conducted  as  follows:  the  contrast  of  the 
stimuli,  f(x,  y)  =f°exp{Acos  (UJy) }  ,  and,  for  each  spatial  frequency  U) , 

A  was  selected  so  that  the  contrasts  of  the  two  stimulus  were  seen  to  be 
the  same.  From  these  selected  A's,  he  computed  the  describing  functions 
for  several  chosen  standard  stimuli. 

An  analysis  of  the  proposed  model  which  predicts  the  variation 
of  the  high  frequency  response  as  a  function  of  contrast  will  now  be 
reported.  For  the  analysis,  it  was  convenient  to  divide  the  system  into 
two  parts  as  shown  in  figure  6(a)  and  investigate  each  separately  then 
combined. 

Consider  the  two  relevant  system  configurations  as  shown  in 
figure  3(a)  and  (b)  Part  1.  It  should  first  be  noted  that  the  systems  are 
not  equivalent.  For  example,  given  a  system  consisting  of  a  log  operation 
followed  by  a  linear  filter,  another  system  consisting  of  a  filter  followed 
by  a  log  operation  can  be  made  equivalent  only  if  the  second  filter 
function  is  made  to  depend  upon  the  input  signal.  However,  for  the  analysis, 
it  was  convenient  to  transform  the  proposed  configuration  into  the  input 
dependent  form  to  obtain  the  exponent-log  cancellation. 

In  the  analysis  it  was  shown  that  the  response  of  the  nonlinear 
combination  of  the  spatial  filter  -  log  did  indeed  produce  a  variation  of  the 
response  with  the  input  signal  contrast.  To  show  this  fact,  the  spatial 
filter  -  log  system  was  first  converted  to  an  equivalent  log  -  input 
dependent  spatial  filter.  The  response  of  this  system  was  expressed 
analylitically  and  numerically  approximated  to  obtain  the  "describing 
functions."  The  results  for  the  low  pass  spatial  filter  and  log  operation 
are  shown  in  figure  4.  The  combined  characteristic  with  the  high  pass 
filter  is  shown  in  figure  15  and  compares  favorable  with  the  experimental 
results  of  Davidson.  Thus,  the  variation  of  the  describing  functions  at 
high  spatial  frequencies  appears  not  to  be  due  to  experimental  equipment, 
but  a  predictable  response  of  the  visual  system. 

The  consequences  of  these  results  also  apply  to  several  other 
theoretical  and  experimental  studies  of  the  visual  system.  For  example, 
a  widely  used  techniques  in  computer  image  processing  is  multiplicative 
filtering  7]  .  In  experimental  work  with  the  log-bandpass  model  the 
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CONTRACT 


Combined  characteristic  spatial  frequency  response  of  the  proposed  model 


contrast  variation  is  predictable.  For  example,  the  relatively  low 
resolution  of  one  lp/mm  viewed  at  40  cm  corresponds  to  about 
7  cycles /degree  and  is  well  above  the  spatial  frequency  at  which  the 
variation  is  noticable.  Further  the  display  of  a  square  wave  image 
display  at  this  frequency  contains  even  higher  frequencies  due  to  the 
harmonics. 

Also,  note  that  the  experiments  such  as  that  performed  by 
Davidson  aimed  at  developing  a  describing  function  of  the  visual  system 
or  those  conducted  by  Stockham  which  illustrated  Mach  band  cancellation 
may  easily  be  performed  with  the  new  model.  Suppose  that  the 
visual  system  is  represented  by  the  cascade  of  a  simple  integrator,  a 
log  nonlinearity  and  a  differentiation.  Then  to  derive  a  describing 
function  of  the  retinal  neural  network  one  may  present  a  sinsoid 

f(y)  =  A  cos  wy 

which  is  processed  by  the  inverse  of  the  log  response  giving 
gj(y)  =  exp  {  f(y)  }  =  exp  {  A  cos  wy} 

and  by  the  inverse  of  the  integration 

d  ,  , 

g 2(y )  =  ~  gj(y)  =  -A  w  sin(wy)  exp  1  A  cos  wyj 

Note  that  the  image  g^fy)  is  now  weighted  in  amplitude  proportionally 
to  the  frequency  of  the  input.  Thus,  as  the  frequency  of  the  input  is 
increased,  the  amplitude  of  the  processed  signal  is  increased.  A  de¬ 
scribing  function  derived  in  th!  s  manner  could  be  called  a  modulation 
transfer  function  and  would  vary  with  the  contrast  of  the  input  signal  in 
a  linear  manner. 

The  preprocessing  required  for  Mach  band  cancellation  is  also 
easily  derived  for  this  simple  example.  Again  starting  with  the  image 


f(y )  =  A  cos  wy 


the  first  step  would  be  an  integration  to  negate  the  effects  of  the  neural 
connections  to  produce 

A 

g^(y)  =  A  cos  (wy)  dy  =  —  sin  (wy)  +  C 
Next  an  exponentation  is  required 


g2(y)  =  exp  {  C]  exp  {  —  sin  wy) 

Finallv,  a  differentiation  is  required 

g3(y)  =  exp  {  C 1  A  sin(  wy)  exp  {  — -  sin  wy) 

The  Image  g^(y)  should  appear  to  the  observer  without  Mach  bands. 

Temporal  Response  .  The  same  descrepancy  between  the  location 
of  elements  responsible  for  the  temporal  frequency  response  in  the 
currently  accepted  model  and  the  proposed  model  may  be  noted.  In  the 
usual  model  the  elements  consists  of  a  nonlinearity  followed  by  the 
temporal  band-pass  filter.  In  the  proposed  model  the  first  element  is  a 
low  pass  temporal  filter,  then  the  nonlinearity  and  the  high  pass  filter. 

Thei  2  is  again  a  physiological  reason  for  placing  the  low  pass 
filter  before  the  nonlinearity.  The  temporal  low  pass  filter  arises  from 
the  30  msec  time  required  for  light  absorbed  by  the  visual  receptors  to 
initiate  the  photochemical  reaction  which  generates  the  neural  pulse. 
Also,  there  is  again  an  inconsistency  encountered  if  one  assumes  the 
log-bandpass  characteristic. 

If  one  postulates  also  the  existence  of  the  excitational  interaction 
among  nerve  cells  as  Yasuda  and  Hiwatashi  [  10]  ,  then  the  neural  net¬ 
work  has  a  band-pass  characteristic.  However,  the  difficulty  in  this 
assumption  is  that,  in  order  to  explain  the  temporal  characteristics  of 
the  H.  V.  S.  ,  one  has  to  assume  that  the  response  time  delay  of  the 
excitational  interaction  is  ignored  whereas  that  of  the  inhibition  is  con- 
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sidered.  Noticing  that  the  excitational  interaction  and  the  inhibitory 
interaction  produce  the  spatial  low-pass  effect  and  the  high  pass  effect, 
respectively,  the  above  statement  means  that  the  low  pass  effect  does 
not  have  a  time  delay  whereas  the  high  pass  effect  does.  Since  the 
low-pass  filter  is  produced  by  the  optical  aberration  of  the  eye  in  the 
model,  the  proposed  model  can  also  overcome  this  difficulty. 

Combined  Spatio-Temporal  Response.  Robson  [  8.)  experimentally 
measured  spatio-temporal  frequency  responses  of  the  visual  system. 

He  used  a  2.  5°  by  2.  5°  grating  target  in  which  the  luminance  at  right  angles 
to  the  bars  was 

L  =  Lq  C  1  +  m(cos2nVy)(cos  2nft)] 

2 

L  was  kept  constant  at  20  cd/m  and  the  value  of  m  at  subjective  dis  - 

0 

appearance  of  bars  was  measured  for  different  frequencies  V  and  f. 

The  inverse  of  that  value  of  m  was  defined  as  the  contrast  sensitivity. 

Since  Robson  did  not  vary  the  input  contrast,  the  predictable  vari¬ 
ation  of  the  resultant  curves  at  high  frequencies  was  of  course  not  ob¬ 
served.  However,  another  interesting  fact  of  the  proposed  model  may 
be  correlated  with  Robson's  results.  This  is  the  separability  of  the 
spatio-temporal  response.  As  pointed  out  by  Budrikis  [  1  J  ,  if  the  human 
visual  system  response  were  product  separable  then  the  family  of  re¬ 
sponse  curves  of  spatial  or  temporal  frequency  with  one  of  the  frequencies 
fixed  would  differ  from  each  other  by  constant  factors  or  in  the  logarith¬ 
mic  plots,  by  constant  displacements  in  the  vertical  direction.  The  ex¬ 
perimental  curves  do  exhibit  this  property  except  at  frequencies  below 
about  5  cycle  j/degree  and  6  cycles/ second.  Note  that  this  is  not  pre¬ 
dicted  by  the  log-bandpass  type  model. 

It  can  be  shown  that  this  effect  is  predictable  from  the 
proposed  model.  The  basic  mechanism  is  that  if  the  output  of  the  non¬ 
linearity  is  at  sufficiently  high  frequencies  so  as  to  be  within  the  passband 
of  the  high  pass  spatio-temporal  filter,  then  the  effect  of  this  high  pass 
filter  is  neglibible.  In  this  situation  the  spatio-temporal  response  is 
product  separable  since  the  low  pass  spatial  and  temporal  filters  are 


independent.  However,  if  the  output  of  the  nonlinearity  contains  low 
frequency  signals  which  are  affected  by  the  stop  band  of  the  high  pass 
filters  then  the  overall  response  is  not  product  separable  since  the 
high  pass  spatio-temporal  filters  are  not  independent. 

Application  to  Color  Vision  .  One  possible  application  to  color 
vision  should  also  be  mentioned.  In  previous  work,  a  difficulty  has 
been  encountered  in  the  study  of  color  balance  which  has  been  attributed 
to  different  shape  nonlinearities  in  the  three  color  receptor  channels. 
However,  it  is  interesting  to  note  that  a  similar  effect  is  predictable 
from  the  proposal  model.  It  is  well  known  that  the  lens  of  the  eye 
produces  chromatic  abberation,  i.e.  different  focal  lengths  are  required 
for  different  spectral  signals.  If  a  color  image  which  may  be  considered 
as  the  sum  of  red,  green  and  blue  component  images,  is  observed,  each 
of  the  component  images  is  seen  blurred  to  a  different  extent.  Thus, 
the  low  pass  filters  in  the  red,  green  and  blue  channels  in  the  proposed 
model,  have  different  cutoff  frequencies.  Even  if  the  nonlinearities  in 
the  three  channels  were  identical,  (and  it  appears  that  there  is  no 
psysiological  evidence  to  support  the  converse)  the  overall  responses  in 
each  channel  would  appear  to  have  different  nonlinearities  [  12]  .  Thus 
it  appears  that  a  correction  for  the  chromatic  abberation  is  essential 
for  quantitative  study  of  color  phenomena. 
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7 •  Image  Processing  Hardware  and  Software  Projects. 


The  image  processing  hardware  and  software  projects  are  devel¬ 
opmental  projects  supportive  of  the  image  processing  research. 

The  hardware  projects  described  in  the  first  report  summarize 
progress  on  the  development  of  real  time  color  image  displays  for  the 
ARPANET.  Also  included  is  progress  on  construction  of  the  real  time 
color  television  recorder  and  playback  system. 

The  software  projects  report  explains  the  image  processing  soft¬ 
ware  system  being  implemented  for  access  over  the  ARPANET. 

7.  1  Hardware  Projects 
Toyone  Mayeda 

A  digital  color  television  system  for  display  of  monochrome  and 
color  television  pictures  received  over  the  ARPANET  has  been  implemented 
ana  in  operation  since  April,  1974.  Development  is  proceeding  on  a  second 
unit  whose  specifications  are  described  below: 

1*  Receive  asynchronous  digital  picture  information  from  the 
ARPANET  TIP  with  brightness  resolution  up  to  256  levels 
and  at  input  rates  up  to  19.  2K  baud. 

2.  Include  a  function  memory  which  can  be  used  to  translate  the 
8  bit  data  words  (from  the  refresh  memory)  with  any  desired 
transfer  curve.  The  function  memory  can  be  remotely  pro¬ 
grammed  from  the  TIP  or  by  local  switch  control. 

3.  Display  a  256  x  256  eight  bit  image,  or  a  six  or  seven  bit  im- 
age  and  a  one  bit  graphic  overlay. 

4.  Use  an  alphanumeric  keyboard  when  available  to  communicate 
with  the  ARPANET  TIP  and  also  to  generate  alphanumeric 
characters  on  the  display  monitor. 

5.  Output  the  monochrome  video  data  and  alphanumeric  charac¬ 
ters  in  composite  RF  format  so  that  it  can  be  displayed  on  an/ 
TV  receiver  using  its  antenna  input. 

The  other  major  hardware  project  under  development  is  a  digital 
magnetic  tape  recorder/playback  unit  for  transfer  of  real  time  color  tele¬ 
vision  signals  to  and  from  a  PDP-10  computer.  The  tape  unit  records  the 
real  time  imagery  in  a  600  ips  mode  and  plays  the  data  back  at  a  1  and  7/8 
ips  rate  which  is  compatible  with  the  PDP-10  channel  capacity.  The  inverse 
process  can  be  performed  to  produce  real  time  television  signals  from 


coded  computer  records.  Presently,  all  control  and  signal  processing  cir¬ 
cuitry  has  been  designed  and  the  high  speed  A/D  has  been  delivered.  De¬ 
livery  of  the  Orion  high  speed  magnetic  tape  recorder  is  expected  in  Decem¬ 
ber,  1974. 

7.2  Software  Projects 
James  M.  Pepin 

The  programming  group  has  worked  on  several  different  projects  in 
the  last  reporting  period.  These  projects  include  support  of  image  proces¬ 
sing  display  devices,  network  programming  and  development  of  an  image 
processing  capability  under  the  Tenex  operating  system. 

The  area  of  progress  in  the  device  support  functions  of  the  Image 
Processing  Lab  include  work  on  support  for  the  microdensitometer  and  im¬ 
plementing  device  routines  for  a  PDP  11/10  acquired  by  the  IPL.  The  sup¬ 
port  for  the  microdensit'i meter  has  been  the  major  area  of  effort  as  this 
device  is  the  most  complex,  from  a  programming  point  of  view,  that  the 
IPL  has  ever  acquired.  The  software  has  progressed  to  the  point  where  it 
is  possible  to  digitize  and  display  data.  The  effort  also  required  that  the 
programming  group  implement  diagnostic  software  to  enable  the  Optronics 
engineers  to  debug  the  device  and  the  IPL  operations  group  to  determine 
the  device's  status.  The  Optronics  is  presently  controlled  by  the  HP2100 
in  the  IPL.  The  plan  is  to  attach  the  device  to  a  PDP  11/10  to  enable  the 
2100  to  perform  other  display  tasks  while  the  Optronics  is  in  use.  Another 
project  is  the  development  of  a  front-end  system  to  the  2100  using  PDP 
11/10's  to  run  the  devices  in  IPL.  This  task  was  taken  to  allow  more  than 
one  device  to  be  in  use  simultaneously.  Previously,  because  of  device 
handshaking  requirements,  the  2100  had  been  only  able  to  run  one  display 
device.  With  the  addition  of  the  11/10's  it  will  be  possible  to  utilize  many 
of  the  devices  in  the  IPL  simultaneously.  The  first  11/10  system  has  been 
implemented  and  a  second  is  on  order.  This  first  system  has  been  inte¬ 
grated  into  the  IPL  operations  and  is  functioning  in  a  production  mode. 

The  largest  area  cf  effort  has  been  the  conversion  of  our  computer 
usage  to  a  KI-10  that  the  ECL  has  acquired.  At  the  present  time  the  KI-10 
has  just  been  installed  and  the  checkout  procedure  of  the  code  developed  is 
just  beginning.  A  great  deal  of  time  has  been  spent  designing  the  inter¬ 
faces  to  the  HP2100  and  related  computers  into  the  IPL.  The  10  is  ex¬ 
pected  to  run  many  of  the  display  devices  with  the  users  interactively.  It 
is  foreseen  that  a  user  can  have  a  Fortran  program  running  on  the  10  in 
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real  time  interacting  with  any  device  in  the  IPL.  This  will  be  implemen¬ 
ted  utilizing  a  PDP  11/40  that  will  be  connected  by  a  DL10  direct  memory 
interface  to  the  KI-10's  memory.  The  KI-10  will  be  running  the  Tenex  op¬ 
erating  system  and  much  time  has  been  spent  familiarizing  the  staff  with 
its  operations  and  usage.  Most  of  the  image  processing  programs  that 
were  implemented  for  the  IBM  360/44  have  been  converted  and  tested  us¬ 
ing  the  PDP-10's  at  ISI. 

The  last  area  that  this  report  will  discuss  is  the  Front  End  nage 
Processing  System  (FEIS).  This  concept  is  to  use  the  KI-10  as  a  front- 
end  to  the  large  "number  crunchers"  on  the  net.  This  concept  has  been 
described  in  the  previous  reports  and  development  of  this  has  entering 
th'‘  early  coding  stages.  The  staff  is  presently  developing  a  set  of  pro¬ 
grams  of  the  IBM  360/91  at  UCLA-CCN.  This  is  the  first  set  of  programs 
being  developed  for  FiES  and  it  is  expected  that  others  will  follow  quickly. 
The  work  on  the  language  interpreter  will  begin  shortly  on  the  KI-10.  Also 
with  the  addition  of  the  KI-10  the  network  handling  functions  will  be  designed 
and  implemented  soon. 
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8.  Image  Processing  Institute  Facilities 

During  the  past  year  the  physical  facilities  of  the  USC  Image 
Processing  Institute  have  been  expanded  considerably.  The  following 
sections  contain  a  summary  description  of  the  present  physical  facili¬ 
ties. 

The  Image  Processing  Institute  is  located  in  Powell  Hall  of 
Information  Sciences,  completed  in  Summer  1973.  Powell  Hall  is  a 
six  story  building  providing  over  100  offices  and  laboratory  facilities. 

The  Institute  laboratories  shown  in  figure  1  include  individual  rooms  for 
digital  image  processing  equipment,  image  display  and  acquisition  de¬ 
vices,  coherent  optics  apparatus,  and  photographic  processors. 

8.1  Image  Processing  Laboratory 

The  Image  Processing  Laboratory  consists  of  four  interconnect¬ 
ing  rooms  housing  image  digitization  and  display  equipment,  a  mini¬ 
computer  system  including  terminals,  magnetic  tape  units,  and  disks 
for  control  of  the  imagp  processing  devices,  a  well  equipped  darkroom, 
and  an  optical  processing  and  holography  laboratory.  The  Image  Pro¬ 
cessing  Laboratory  computer  system  is  interconnected  with  the  Biomed¬ 
ical  Image  Processing  Laboratory  and  the  image  processing  computer 
system  located  in  the  Engineering  Computer  Laboratory. 

Figure  2  contains  a  block  diagram  of  the  Image  Processing  Labora¬ 
tory  computer /controller  system.  This  system  consists  of  a  Hewlett- 
Packard  Model  2100  digital  computer  with  a  variety  of  computer  peri¬ 
pherals  accessed  through  an  input-output  controller.  The  I/O  controller 
also  acts  as  the  communciations  link  with  the  image  processing  com¬ 
puter  system  of  the  Engineering  Computer  Laboratory  and  the  Biomedical 
Image  Processing  Laboratory. 

The  Laboratory  image  display  and  digitization  devices  are  listed 
in  figure  3  and  shown  in  the  photographs  of  figures  4  to  6.  With  these 
devices  it  is  possible  to  digitize  monochrome  and  color  transparencies 
ranging  in  size  from  16mm  up  to  full  size  X-rays  and  pr:nts  of  up  to 
8  in.  x  10  in.  in  size.  Hard  copy  monochrome  and  color  transparencies 
and  prints  can  be  produced  in  a  wide  range  of  sizes  up  to  8  in.  x  10  in. 
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Figure  8-1  Floor  plan  of  powell  hall  laboratories 
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Figure  8-2  Image  processing  laboratory  computer/controller  facilities. 


Figure  8-3  lmage  processing  digitization  and 
display  devices. 


a)  Digital  Process  Room 


b)  Left  Hand  View  of  Im 
Equipment  Room 


c)  Richt  Hand  View  of  lmace 
Equipment  Room 


Figure  8-4.  Facilities  of  the  Digital  Process  Room  and  Image  Equipment  Roan 


a)  Muirhead  Simultaneous 
Color  Scanner 


b)  IER  Flying  Spot 
Scanner 


c)  Dicomed  14"  x  17' 
Scanner 


Figure  8-5.  Some  Devices  Used  for  Computer-Image  I/O. 


a)  Real  Time  Display 


b)  Flatbed  Microdensitometer 


c)  USC  Developed  Display 


Figure  8-6.  Some  Devices  Used  For  Computer-Image  I/O. 


I 


The  IER  flying  spot  scanner  is  used  to  obtain  Polaroid  prints 
and  transparencies  of  images.  Resolution  is  up  to  1024  x  1024  pixels 
with  8  bits  of  intensity  quantization.  Scan  times  vary  from  10  to  60 
sec.  depending  upon  resolution.  Color  photographs  are  obtained  by 
sequentially  scanning  the  white  CRT  phosphor  through  color  filters. 


The  Muirhead  color  facsimile  device  is  a  two  unit  drum  trans¬ 


mitter  and  receiver.  Its  capabilities  include  100  line/inch  resolution 
with  simultaneous  8  bit/primary  quantization.  The  scan  time  is  12  min. 
for  an  8  in.  x  10  in.  print. 

The  Dicomed  14  in.  x  17  in.  X-ray  and  large  transparency 
scanner  is  an  image  dissector  camera  mounted  in  a  light  box.  It 
is  capabi’  of  digitizing  with  8  bits  of  accuracy  and  a  spatial  resolution 
of  up  to  2048  pixels. 

The  Optronics  precision  flat  bed  scanning  microdensitometer 
scans  and  records  16,  35,  and  70mm  color  or  monochrome  roll  film 


on  a  registered  transport  under  computer  control  as  well  as  prints  or 


transparencies  up  to  10  inches  square.  Spatial  accuracy  is  _+  1  micron 
with  apertures  as  small  as  10  microns.  Photographic  density  or  trans¬ 
mittance  is  digitized  with  12  bits  and  the  display  is  driven  by  10  bits. 
The  scanning  velocity  is  8  inches/sec  and  the  unit  is  capable  of  digit¬ 
ization  over  a  0  to  4  specular  density  range.  Color  digitization  and 
display  is  performed  sequentially  under  computer  control. 

Real  time  display  of  monochrome  and  color  digital  images  is 
accomplished  using  an  Aerojet  General  display  device.  This  device 


utilizes  a  standard  shadow  mask  CRT  with  576  horizontal  and  525  vertical 


lines  of  resolution.  A  digital  disk  is  included  with  the  device  which 
makes  possible  a  refresh  rate  of  60  field s  /  second  at  64  quantization 
levels  for  each  of  the  red,  green,  and  blue  primaries. 

A  real  time  display  device  designed  to  be  plugged  directly  into 
the  ARPA-TIPhas  been  developed  by  the  USC  Image  Processing  Institute. 
A  keyboard  terminal  controls  the  selection  of  imagery  from  the  net  *o  be 
displayed  on  the  device. 
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Figure  7  contains  a  block  diagram  of  a  real  time  digital  color 
television  recording  and  playback  system  being  developed  by  the  USC 
Image  Processing  Institute.  In  this  system  the  color  signal  output  of 
a  television  camera  or  video  tape  recorder  is  sampled  and  digitized 
by  three  analog-to-digital  converters  at  a  rate  compatible  with  U.  S. 
television  standards.  The  digitized  signal  is  then  recorded  on  the 
high  speed  Orion  Titan  digital  tape  recorder  shown  in  figure  8.  This 
tape  recorder  is  capable  of  recording  and  playing  back  at  rates  of  up  to 
10  Mbs  on  each  of  28  data  tracks  for  a  total  of  280  Mbs  storage  rate.  A 
single  reel  of  tape  can  record  over  one  minute  of  real  time  television 
data  at  this  rate.  The  recorder  is  also  capable  of  playing  back  at  a  320:1 
speed  reduction  to  match  the  channel  rate  of  the  PDP-11  minicomputer 
interface  to  the  PDP-10.  In  the  output  mode  the  PDP-10  will  output  images 
for  recording  on  the  Orion  at  the  slow  rate;  and  the  Orion  will  then  play¬ 
back  at  its  higher  speed  to  reproduce  a  reconstructed  real  time  tele¬ 
vision  signal. 

8.2  Engineering  Computer  Laboratory 

The  Engineering  Computer  Laboratory  is  based  upon  a  dual  pro¬ 
cessor  computer  system  shown  in  figure  9  consisting  of  Digital  Equipment 
Corp.  PDP  Model  KI10  processor  for  time  share  computing  and  an  IBM 
360/44  processor  for  moderate  size  batch  computing.  Both  machines  are 
host  computers  on  the  ARPA  computer  network.  The  PDP-10  utilizes 
the  TENEX  operating  system  while  the  360/44  runs  under  a  USC  developed 
operating  system  that  incorporates  the  VICAR  image  processing  software 
system  developed  at  the  Jet  Propulsion  Laboratory.  The  PDP-10  and 
360/44  are  connected  together  and  to  the  ARPANET-TIP  through  a 

communications  controller.  The  controller  provides  character 
manipulation  for  ARPANET  transmission  and  interfacing  to  the  University 
Computer  Center  IBM  3  70/158,  the  Image  Processing  Laboratory  image 
digitization  and  display  devices  including  the  real  time  color  TV  digi¬ 
tizer  and  recorder,  and  to  the  Adage  AGT-10  computer  graphics  system. 

The  Image  Processing  Institute  is  developing  a  front  end  image 
processing  system  (FEIS)  on  the  PDP-10  which  will  be  available  to  any 
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image  processing  user  at  USC  or  over  the  ARPANET.  Specifically 
FEIS  will  make  use  of  the  existing  ARPANET  structure  on  a  dynamic 
basis  to  allocate  the  image  scanning,  number  crunching,  dispiay,  trans¬ 
porting,  etc.  resources  of  the  network  for  specific  image  processing  job 
requirements.  The  system  will  be  transparent  to  the  user  in  the  sense 
that  the  user  need  not  possess  intimate  knowledge  of  resources  available 
on  the  network,  nor  will  the  user  be  required  to  have  an  intimate  know¬ 
ledge  of  existing  image  processing  software,  protocols,  or  network 
procedures.  Figure  10  indicates  a  brief  overview  of  how  such  a  system 
might  operate.  From  the  figure  it  is  apparent  that  three  types  of  infor¬ 
mation  flow  are  under  control  of  FEIS.  Naturally  there  would  be  TELNET 
(network  communication)  and  RJE  (remote  job  entry)  command  and  control 
instructions  in  which  FEIS  would  command  large  image  processing  jobs 
to  be  run  on  the  number  crunchers  of  the  network.  However  a  typical 
user  may  execute  his  job  in  ANIPIL  (ARPANET  Image  Processing 
Interpretive  Language)  directly  to  the  USC  front  end  system.  Finally 
image  flow  would  be  routed  through  the  most  efficient  path  possible  to 
preserve  the  precious  bandwidth  needed  for  large  image  arrays  and  to 
establish  and  minimize  an  archival  image  storage  facility  for  user's 
future  reference. 

8.3  Optical  Information  Processing  Laboratory 

As  a  compliment  to  research  efforts  in  digital  image  processing 
now  being  conducted  at  the  USC  Image  Processing  Institute,  an  Ootical 
Information  Processing  Laboratory  (OIPL)  has  been  established  within 
the  Institute  and  the  Department  of  Electrical  Engineering.  The  equip¬ 
ment  available  in  this  laboratory  combined  with  facilities  for  photography 
and  for  digital  scanning,  computing  and  display  which  are  already  in 
operation  provide  a  unique  capability  for  research  in  the  fields  of  optical 
and  hybrid  optical /digital  information  processing. 

The  OIPL  is  physically  located  in  two  rooms  on  the  first  floor 
of  Powell  Hall,  adjacent  to  the  digital  image  processing  facilities.  One 
room  contains  a  steel  topped  honeycomb  optical  stable  table,  4'  x  10'  x  8" 
in  size,  along  with  a  computer-controlled  film  transport  and  ROSA 
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ARPA  picture  -  typical  mode  of  operation. 


(recording  optical  spectrum  analyzer)  photodetector  and  data  collection 
system.  This  optical /digital  pattern  recognition  system  also  includes 
a  15  mW,  He-Ne  laser  source,  mirrors,  collinating  optics,  pinhole 
spatial  filter,  positioners  and  mounts  which  can  be  rearranged  into  dif¬ 
ferent  system  configurations. 

Another  room,  located  at  the  rear  of  the  IPI  computer  room,  is 
entirely  devoted  to  general  research  in  optical /digital  processing.  The 
room  is  light-tight,  and  different  types  of  experimental  arrangements 
are  possible.  Three  optical  tables  are  available;  one  is  a  4'  x  12’  x  12" 
steel-honeycomb  table;  the  others  are  granite  surface  plate  tables 
4'  x  4'  x  8"  and  2'  x  8'  x  12".  Both  granite  tables  have  innertube  inflat¬ 
able  suspension  systems  for  vibration  isolation  in  interferometric  systems. 
The  lab  has  both  15  mW  and  a  50  mW  He-He  coherent  laser  sources, 
i  m  and  3  m  optical  rails  and  mounts,  and  a  limited  assortment  of 
positioners,  mirrors,  lenses,  flats,  prisms,  rulings,  collinators  and 
other  optical  equipment  for  experimentation.  Measuring  equipment  in¬ 
cludes  transmission  and  reflection  densitometers,  a  radiometer /photo¬ 
meter  and  an  electrical  shutter  and  timer.  Photographic  facilities  for 
processing  the  high  resolution  plates  and  films  generally  used  in  optical 
systems  are  now  being  developed. 

All  these  facilities  can  be  utilized  and  adapted  for  research  in 
real-time  optical /digital  systems  making  use  of  fast  controllable  deflectors, 
modulators,  and  detectors.  Several  projects  low  underway  make  use  of 
optical  masks  or  screen  produced  digitally  on  IPI  displays.  Although 
work  in  these  areas  is  still  at  an  early  stage,  research  at  these  facili¬ 
ties  should  develop  the  potential  of  optical,  digital,  and  hybrid  information 
proce  s  sing. 

8.4  Biomedical  Image  Processing  Laboratory 

The  Biomedical  Image  Processing  Laboratory  has  been  establish¬ 
ed  in  augmentations  of  the  Image  Processing  Laboratory  to  provide  re¬ 
search  facilities  in  this  specific  area  of  research.  Presently,  the  Labora¬ 
tory  contains  an  Adage  Model  AGT-10  computer  graphics  system  and 
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several  computer  terminals  connected  to  the  TIP.  The  Adage  possesses 
stand-alone  graphics  capability  for  interactive  graphics  computing,  and 
is  also  connected  to  the  IBM  360/44  computer  for  operations  requiring 
greater  computer  power. 

A  Recognition  Systems  Incorporated  computer/laser  optical 
processing  system  shown  in  figure  11  is  installed  in  the  Laboratory.  This 
unit  contains  a  coherent  imaging  apparatus  that  produces  a  Fraunhofer 
diffraction  light  pattern  of  a  photographic  transparency.  The  diffraction 
pattern  is  sensed  by  an  optical  detector  containing  64  elements  that 
sample  spatial  frequency  and  amplitude  components  of  the  light  pattern. 
The  64  outputs  are  then  processed  by  analog  circuitry  and  fed  to  the 
HI  -2100  of  the  Image  Processing  Laboratory  for  real  time  data  reduction. 
Similar  equipment  has  been  used  quite  successfully  by  USC  staff  members 
for  photographic  texture  detection  in  the  development  of  automatic  equip¬ 
ment  for  the  diagnosis  of  coal  miner's  pneumoconiosis  and  is  directly 
applicable  to  reconnaissance  and  other  DoD  related  imagery. 
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